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Abstract  We  develop  two  characteristic  methods  for  the  solution  of  the  linear  advection  dif¬ 
fusion  equations  which  use  a  second  order  Runge-Kutta  approximation  of  the  characteristics 
within  the  framework  of  the  Eulerian-Lagrangian  localized  adjoint  method.  These  methods  nat¬ 
urally  incorporate  all  three  types  of  boundary  conditions  in  their  formulations,  are  fully  mass 
conservative,  and  generate  regularly  structured  systems  which  arc  symmetric  and  positive  def¬ 
inite  for  most  combinations  of  the  boundary  conditions.  Extensive  numerical  experiments  arc 
presented  which  compare  the  performance  of  these  two  Runge-Kutta  methods  to  many  other  well 
perceived  and  widely  used  methods  which  include  many  Galerkin  methods  and,  high  resolution 
methods  from  fluid  dynamics. 
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1  Introduction 

Advection-diffusion  equations  are  an  important  class  of  partial  differential  equations  that 
arise  in  many  scientific  fields  including  fluid  mechanics,  gas  dynamics,  and  atmospheric 
modeling.  These  equations  model  physical  phenomenon  characterized  by  a  moving  front. 
In  fluid  dynamics,  for  example,  the  movement  of  a  solute  in  ground  water  is  described  by 
such  an  equation.  Since  these  equations  normally  have  no  closed  form  analytical  solutions, 
it  is  very  important  to  have  accurate  numerical  approximations.  When  diffusion  dominates 
the  physical  process,  standard  finite  difference  methods  (FDM)  and  finite  element  methods 
(FEM)  work  well  in  solving  these  equations.  However,  when  advection  is  the  dominant  pro¬ 
cess,  these  equations  present  many  numerical  difficulties.  In  fact,  standard  finite  element  and 
finite  difference  methods  produce  solutions  which  exhibit  non-physical  oscillations,  excessive 
numerical  diffusion  which  smears  out  sharp  fronts,  or  a  combination  of  both  [21,  35]. 

Many  specialized  schemes  have  been  developed  to  overcome  the  difficulties  mentioned. 
One  class  of  these  methods,  usually  referred  to  as  the  class  of  Eulerian  methods ,  uses  an 
Eulerian  fixed  grid  and  improved  techniques,  such  as  upstream  weighting ,  to  generate  more 
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accurate  approximations.  Most  of  the  methods  in  this  class  are  characterized  by  ease  of  for¬ 
mulation  and  implementation,  however  their  solutions  suffer  from  excessive  time  truncation 
errors.  Moreover,  they  put  a  strong  limitation  on  the  Courant  number  and  hence  require  very 
small  time  steps  to  generate  stable  solutions.  Among  these  methods  are  the  Petrov-Galerkin 
FEM  methods  [1,  4,  9,  55]  which  are  improvements  over  the  standard  Galerkin  FEM  that 
incorporate  unwinding  in  the  space  of  the  test  functions.  Also  included  in  the  class  of  Eule- 
rian  methods  are  the  streamline  diffusion  FEM  methods  (SDM)  [5,  19,  26,  31,  32,  33,  35,  36] 
and  the  continuous  and  discontinuous  Galerkin  methods  (CGM,  DGM)  [23,  34,  37,  42,  43]. 
The  SDM  improve  over  the  standard  space-time  Galerkin  FEM  by  adding  a  multiple  of 
the  (linearized)  hyperbolic  operator  of  the  problem  considered  to  the  standard  test  func¬ 
tions.  Thus  they  add  numerical  diffusion  only  in  the  direction  of  the  streamlines.  The  SDM 
formulations  have  a  free  parameter  which  determines  the  amount  of  diffusion  applied  and 
therefore  has  a  great  effect  on  the  accuracy  of  these  methods.  In  practice,  this  parameter 
should  be  large  enough  to  avoid  oscillations  in  the  solution,  but  not  too  large  to  damp  the 
solution.  A  clear  choice  of  this  parameter  is  not  known,  in  general,  and  is  heavily  problem 
dependent.  The  CGM  and  DGM  are  well  suited  for  purely  hyperbolic  equations  and  recently 
have  been  extended  to  solve  the  advection-diffusion  equations.  They  are  space-time  explicit 
methods  in  which,  starting  with  the  initial  time  solution  and  Dirichlet  data  at  the  inflow 
boundary,  one  would  successively  iterate  over  the  elements  of  a  quasi-uniform  triangulation 
of  the  space-time  domain  in  an  order  consistent  with  the  domain,  solving  a  local  system 
over  each  element.  In  addition  to  the  methods  mentioned  above,  the  class  of  Eulerian  meth¬ 
ods  includes  the  high  resolution  methods  in  fluid  dynamics  such  as  the  Godunov  methods, 
the  total  variation  diminishing  methods  (TVD),  and  the  essentially  non-oscillatory  methods 
(ENO)  [10,  11,  12,  18,  24,  27,  46,  47,  48,  49].  These  methods,  as  well  as  the  CGM  and 
DGM,  are  well  suited  for  advection-diffusion  equations  with  small  diffusion  coefficients  and 
in  general  impose  an  extra  stability  restriction  on  the  size  of  the  time  step  taken  based  on  the 
magnitude  of  this  coefficient.  Therefore,  they  are  very  sensitive  to  changes  in  the  diffusion 
coefficient,  which  in  practical  problems  is  likely  to  have  large  values  at  certain  points. 

Another  class  of  methods,  usually  known  as  characteristic  methods ,  makes  use  of  the 
hyperbolic  nature  of  the  governing  equation.  These  methods  use  a  combination  of  Eulerian 
fixed  grids  to  treat  the  diffusive  component,  and  Lagrangian  coordinates  by  tracking  particles 
along  the  characteristics  to  treat  the  advective  component.  Included  in  this  class  are  the 
Eulerian-Lagrangian  methods  (ELM),  the  modified  methods  of  characteristics  (MMOC), 
and  the  operator  splitting  methods  [2,  13,  17,  20,  21,  25,  30,  38,  39,  40,  41,  50,  56].  These 
methods  have  the  desirable  advantage  of  alleviating  the  restrictions  on  the  Courant  number, 
thus  allowing  for  large  time  steps.  Furthermore,  the  Lagrangian  treatment  in  these  methods 
greatly  reduces  the  time  truncation  errors  which  are  present  in  Eulerian  methods.  On  the 
other  hand  these  methods  have  difficulty  in  conserving  mass  and  treating  general  boundary 
conditions.  The  Eulerian-Lagrangian  localized  adjoint  methods  (ELLAM)  were  developed 
as  an  improved  extension  of  the  characteristic  methods  that  maintains  their  advantages  but 
enhances  their  performance  by  conserving  mass  and  treating  general  boundary  conditions 
naturally  in  their  formulations.  The  first  ELLAM  formulations  were  developed  for  constant 
coefficient  advection-diffusion  equations  [7,  44],  The  strong  potential  that  these  FEM  based 
formulations  and  their  numerical  results  have  demonstrated  have  led  to  the  development 
of  additional  formulations  for  variable-coefficient  advection  diffusion  equations  [45,  51]  and 
for  nonlinear  equations  [14]  as  well  as  finite  volume  formulations  [8,  28].  However,  because 
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the  characteristics  for  variable-coefficient  advection-diffusion  equations  cannot  be  tracked 
exactly  in  general,  many  characteristic  and  ELLAM  methods  that  have  been  developed  use 
a  backward  Euler  approximation  for  these  characteristics  due  to  its  simplicity  and  stability. 
These  (backward  Euler)  schemes  are  second  order  accurate  in  space  but  only  first  order 
accurate  in  time. 

An  ELLAM  based  formulation  which  uses  a  second  order  approximation  for  the  char¬ 
acteristics  was  developed  recently  for  first-order  advection-reaction  equations  [54].  This 
formulation  was  shown  to  be  of  second  order  in  both  space  and  time.  Unlike  the  first-order 
equations,  similar  treatment  for  the  advection-diffusion  equations  is  more  problematic  since 
one  needs  to  treat  the  diffusive  component  and  its  partial  derivatives  (with  respect  to  the 
rectangular  and  characteristic  coordinates)  carefully  along  the  characteristics  to  produce 
systems  having  desirable  structure.  In  addition,  boundary  treatment  is  more  involved  be¬ 
cause  at  both  the  inflow  and  outflow  boundary  data  can  be  specified  in  many  forms  which 
then  need  to  be  incorporated  into  the  formulations.  In  this  paper  we  develop  characteris¬ 
tic  methods  for  the  one-dimensional  linear  advection-diffusion  equations,  based  on  a  second 
order  Runge-Kutta  approximation  of  the  characteristics.  We  present  a  backward  tracking 
(BRKC)  and  a  forward  tracking  (FRKC)  Runge-Kutta  characteristic  scheme,  both  of  which 
are  mass  conservative  and  incorporate  boundary  conditions  naturally  in  their  formulations. 
These  methods,  which  can  be  thought  of  as  generalizations  of  the  ELLAM  schemes,  generate 
tridiagonal  (regularly  structured  in  multi-dimensions)  matrices  which  are  symmetric  (except 
at  the  inflow  boundary)  and  positive  definite  that  can  be  solved  efficiently.  Moreover,  we 
provide  the  results  of  some  extensive  numerical  experiments  which  compare  the  performance 
of  the  two  methods  developed  to  many  of  the  methods  mentioned  above. 

This  paper  is  organized  as  follows.  In  Section  2  we  develop  a  reference  equation  based  on 
exact  characteristic  tracking.  In  Section  3  we  present  the  two  characteristic  schemes  BRKC 
and  FRKC  which  are  based  on  a  backward  tracking  and  a  forward  tracking  algorithm, 
respectively.  We  also  give  a  detailed  description  of  boundary  treatment,  in  addition  to  some 
numerical  experiments  which  demonstrate  the  order  of  convergence  of  the  two  schemes.  In 
Section  4  we  give  a  brief  description  of  some  well  studied  and  widely  used  methods  which  are 
known  to  give  good  approximations  to  the  advection  diffusion  equations.  Section  5  contains 
the  results  of  the  numerical  experiments  that  compare  the  performance  of  the  two  schemes 
developed  in  Section  3  and  the  other  methods  described  in  Section  4. 

2  Development  of  the  Method 

We  consider  the  one-dimensional  linear  variable-coefficient  advection-diffusion  equation 

Cu  =  ut  +  (V(x,t)u  -  D(x,t)ua,)x  -  f(x.l).  x  G  (a,  b),  i  G  [0,T],  . 

u(x,  0)  =  ua(x),  rr  G  [a,  6], 

where  V(x,t)  and  D(x,t )  are  the  velocity  field  and  the  diffusion  coefficient,  respectively. 
D(x,t )  is  assumed  to  be  positive  throughout  the  domain.  To  simplify  our  presentation, 
we  also  assume  V(a,t )  and  V(b,t )  are  positive,  i.e.  we  set  x  =  a  and  x  =  b  to  be  the 
inflow  and  outflow  boundaries,  respectively.  In  many  advection-dominated  applications, 
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I  D(x,  f)|  <<  | V" (rc, t) | ,  therefore,  methods  devised  to  solve  equation  (2.1)  should  handle  this 
case  accurately.  We  consider  general  boundary  conditions  with  any  combination  of  Diriclilet, 
Neumann,  or  Robin  conditions  at  the  inflow  and  outflow  boundaries.  The  function  g(t)  is 
used  to  denote  the  time  dependence  to  specify  the  boundary  condition  at  x  =  a  and  h(t) 
for  the  outflow  boundary  at  x  =  b.  The  Dirichlet,  Neumann,  and  Robin  conditions  are 
formulated  by  the  requirements 


u(c,  t ) 

=  0i  (*), 

t  e  (o, T], 

( Dux)(c ,  t) 

=  02  (t), 

1  c  (0.7’], 

(2.2) 

■  Dux)(c ,  t) 

=  03  W, 

t  e  (0, T], 

respectively,  where  <f>  =  g  when  c  =  a  and  <j>  =  h  when  c  =  b. 


2.1  Variational  Formulation  and  Characteristic  Curves 

The  domain  of  problem  (2.1),  is  [a,  b]  x  [0,  T]  which  we  partition  in  space  and  time  as  follows: 


0  =  t°  <  t1  <  ...<  tN  =  T, 

a  =  x'q  <  x "  <  . . .  <  x'j'  =  b,  »  —  ().! . V  . 


(2.3) 


for  positive  integers  N,  In  (n  =  0, 1,  ...,  N).  We  utilize  space-time  test  functions  that  vanish 
outside  [a,  b]  x  (tn,tn+1],  which  enables  us  to  concentrate  on  one  time  period  (tn,.in+1].  Our 
test  functions  are  discontinuous  in  time  and  allow  for  different  meshes  at  different  time  peri¬ 
ods.  For  notational  simplicity,  we  suppress  the  temporal  index  on  the  grid.  The  variational 
formulation  of  equation  (2.1)  on  the  domain  Q  =  [a,  6]  x  (/." , /;"+ 1  ].  obtained  by  multiplying 
that  equation  by  a  test  function  w  (which  we  describe  in  more  detail  below)  and  integrating 
by  parts,  becomes 


/  u(x,f1+1)w(x,tn+1)  dx  +  /  /  Duxwx  dxdt 

Ja  Jtn  Ja 

^in+1  i>tn+1  ,-h 

+  /  (Vu  —  Dux)w  | ba  dt  —  /  u(wt  - f-  Vwx) 

Jt"  Jtn  Ja 

pb  fitn~^1  fib 

=  /  a(xii"  )w(x.  t"_  )  dx  +  /  /  fw  dxdt, 

Ja  Jtn  Ja 


dxdt 


(2.4) 


where  we  use  the  notation  w(x,t+)  =  limf^f..  w (x,t)  due  to  the  discontinuity  of  w(x,t )  at 
time  tn. 


The  principle  of  the  localized  adjoint  methods  (LAM)  requires  the  test  functions  to  be 
chosen  from  the  solution  space  of  the  homogeneous  adjoint  equation 

CT  w  =  —wt  —  V(x,t)wx  +  (D(x,t)w =  0.  (2.5) 

However,  the  solution  space  of  this  equation  is  infinite  dimensional.  Thus  we  need  to  split 
this  equation  to  determine  a  finite  number  of  test  functions.  In  fact,  different  splittings  lead 


4 


to  different  approximation  schemes  (see  [7,  52]).  Due  to  the  Lagrangian  nature  of  the  exact 
solution,  a  natural  splitting  is 


wt  +  V(x,  t)wx  =  0, 
(J)  (  X  ,  I  )  ((',  ) -  0. 


(2.6) 


The  ELLAM  generalizes  the  framework  of  the  localized  adjoint  methods  by  requiring  the 
test  functions  to  satisfy  the  first  equation  in  (2.6)  exactly  or  even  approximately  (so  that 
the  last  term  on  the  left  hand  side  of  equation  (2.4)  vanishes  exactly  or  approximately). 
However,  the  ELLAM  does  not  require  the  second  equation  in  (2.6)  to  be  satisfied,  thus  one 
does  not  have  to  choose  the  test  functions  w(x,t )  to  be  of  a  complicated  form  in  space.  In 
our  formulation,  we  choose  the  test  functions  to  be  piecewise  linear  in  space,  as  in  the  linear 
standard  FEM,  and  to  be  constant  along  characteristics  which  we  discuss  below. 


The  characteristic  curves  of  (2.1)  are  given  by  the  solutions  of  initial  value  problems  for 
the  ordinary  differential  equation 

§  =  v<*‘>-  (2-7) 

We  denote  the  characteristic  curve  emanating  from  a  given  point  (x,t)  with  t  £  [f",f"+1],  by 

V  =  X(0]  x,t),  (2.8) 

where  8  is  the  time  position  parameter  along  that  characteristic.  Furthermore,  we  introduce 
below  some  notation  that  is  described  by  the  following  relations 


X 

=  X(tn+1]x, 

X* 

=  X(t":  x.  1." 

b'lf) 

=  X(tn]b,t), 

a 

=  A’ (  f  (.r):  x 

b 

=  X(t(x)]x, 

(2.9) 


We  define  x  and  x*  as  the  head  (using  forward  tracking)  and  the  foot  (using  backward 
tracking),  respectively,  of  the  characteristics.  The  foot  of  a  characteristic  with  head  on  the 
outflow  boundary  is  denoted  by  b*(t)  for  t  £  [tn,tn+1].  We  also  define  t(x)  and  t*(x)  as 
the  exit  times  of  the  characteristics  X(8]x,tn)  and  X(8]  x,tn+1)  at  the  outflow  and  inflow 
boundaries,  respectively  (see  Figure  l.b).  The  general  time  increment  over  the  domain  may 
then  be  written  as 

A t(x)  :=  t(x)  —  t*(x), 

where  we  define  t*(x)  =  t"  for  characteristics  with  feet  not  on  the  inflow  boundary,  and 
similarly  t(x)  =  t”+1  for  characteristics  with  heads  not  belonging  to  the  outflow  boundary. 
By  implicitly  differentiating  the  fourth  relation  in  (2.9)  for  t*(x)  with  respect  to  x,  and 
the  relation  X(8\b,t )  =  b  +  Jte  V(X(s\b,t),  s)ds,  for  characteristics  X(6\b,t )  with  8  £  [f\t] 
originating  at  points  (b,t)  on  the  outflow  boundary,  with  respect  to  t,  we  get  the  following 
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equations  (for  partial  derivatives  of  the  characteristics) 


Xx(t*(x)]x,tn+1) 

Xt(8-b,t)\0=t 

in  respective  order.  In  the  following  section,  we  consider  approximations  of  the  characteris¬ 
tics,  defined  by  second  order  Runge-Kutta  approximations. 

2.2  Reference  Equation 

After  we  have  introduced  the  notion  of  the  characteristics,  we  are  able  to  find  a  reference 
equation  by  expanding  certain  of  the  integrals  in  the  variational  formulation  (2.4).  First  we 
assume  that  we  can  exactly  track  the  characteristics,  thus  the  same  reference  equation  is 
generated  if  we  use  either  forward  or  backward  tracking  algorithms. 


=  -V(a,t*(x)) 
=  ~ V(b,t ), 


dtx  (x) 
dx 


(2.10) 


2.2.1  Treatment  of  the  Source  Term 


We  start  with  the  source  term  (second  term  on  the  right  hand  side  of  equation  (2.4))  which 
we  break  as  follows, 


f  [  fwdxdt=  [[  fwdyds+  [[  fwdyds+  [[  fwdyds.  (2.11) 

m  Ja  JJch  JJch  JJ  n3 


Here  the  region  fR  represents  all  points  which  he  on  characteristics  with  feet  on  the  inflow 
boundary,  0:>  represents  all  points  on  characteristics  with  heads  on  the  outflow  boundary, 
and  finally  O2  represents  the  remainder  of  Q  (see  Figure  l.a).  For  clarity  of  presentation, 
we  use  the  variable  y  to  represent  the  spatial  coordinate  of  any  point  in  O,  and  reserve  x 
for  points  on  the  spatial  mesh  of  O  at  time  tn  or  tn+1,  representing  either  heads  or  feet  of 
characteristics.  Similarly  we  let  s  denote  the  temporal  variable  instead  of  t,  which  we  reserve 
for  the  temporal  coordinate  of  heads  or  feet  of  characteristics  which  lie  on  either  the  inflow 
or  outflow  boundary.  The  general  time  increment  can  also  be  described  for  points  x  at  time 
tn+1  or  time  t"  in  respective  order  by  A t(I\x  )  and  A ti'°\  x)  which  are  defined  as  follows 


At^  (x)  =  tn+1  —  t*(x),  x£[a,b\, 
A t{°\x)  =  t(x)  —  V\  x  e[a,  6], 


(2.12) 


where  t*(x)  and  t(x)  extend  to  i,n  and  tn+1,  respectively,  if  the  characteristics  that  define 
them  lie  in  02.  With  the  change  of  variable  y  =  y(x,  s )  =  X(s]  x,  t”+1)  for  the  first  integral 
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on  the  right  hand  side  of  (2.11),  we  obtain 


[Jo  f(y,s)  w(y,s)  dyds 

=  f  I  f(X(s]x,t'1+1),s)w(X(s-1x,tn+1),s)Xa,(s]x,tn+1)dsdx 
Ja  Jt*(x ) 

=  J  M  0  ^  (/(M”+1)  w(x,tn+1 )  + 


/(a,  t*(x))  w(a,  t*(x))  Xx(t*(x)\  x,  f"+1)j  dx  +  Ef,Q1(w) 

r  At{T)  ^f(x,tn+i)  w(x,tn+i )  d# 

./a  2 

L 


t."  + 1  (+n+ 1 


+ 


(t"+1  -  f) 


/(a.  1)  w(a,  t)  V (a, t)  dt  +  i?/i£il  (w). 


(2.13) 


where  in  the  second  equality  we  used  a  trapezoidal  approximation  for  the  inner  integral 
whose  error  Ef (w)  is  given  below.  In  the  last  equality  we  moved  the  second  integral  to 
the  inflow  boundary  by  the  change  of  variable  a  =  X(t*(x)  \  x,  tn+1 )  and  the  first  equation  in 
(2.10).  The  trapezoidal  error  term  introduced  is  given  by 


Ef,nt  (w)  = 

E  rtn+1  (t*(x)  -  s)(tn+1  -  s)  d2 

Ja  Jt*(x )  2  d.S 2 


Xx(s;  tn+1)f(X(s]  x,  tn+1),  s)  w(x,tn+1)  dsdx , 


(2.14) 

since  w  is  constant  along  characteristics.  With  the  same  change  of  variable,  y  =  X (s;  x.  t"+1), 
and  similar  treatment  as  in  (2.13),  we  expand  the  second  integral  on  the  right  hand  side  of 
equation  (2.11)  to  get 

//  f(y,s)w(y,s)dyds 
JJa2 

=  f  f  f(X(s]  %i”+1),  s)  w(X(s]  x,  tn+1)T.is)  Xx(s;  x,  tn+1)  dsdx 

Ja  Jtn 

rh  At  (  dx*  \ 

=  I  -y  lf{x,tn+1)  w(x,tn+1 )  +  fix" .  r )  w(x*,tl)  —J  dx  +  Ef,  n2(w) 

=  /  —f(x,  tn+1 )  w(x,  tn+1)  dx  +  — -f(x ,  tn )  w(x.  t’l )  dx  +  /'/.q.,  (w), 

Ja  J  Ja  Z 


(2.15) 


where  we  have  used  the  fact  that  =  Xx(tn\  x, tn+1).  In  this  expression  the  trapezoidal 
error  term  is  given  by 


Ef.n2(w)  - 

\  72 

—  — -  x,  tn+1)f(X(s]  §|,  tn+1),  s)j  w(x,  tn+1)  dsdx. 

(2-16) 

The  treatment  of  the  source  term  over  fi3  is  slightly  different  since  for  the  characteristics, 
we  backtrack  from  points  originating  on  the  outflow  boundary.  With  the  change  of  variable 


/7 


Lfl-fl 


(r  -  s)(tn+1  - 
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y  =  X(s;  b,t ),  the  third  term  on  the  right  hand  side  of  (2.11)  becomes 

//  f(y,s)w(y,s)dyds 

././n3 

rtn+1  rl  ......  ....  ..  . .. ,  .  , 

=  —  /  /  f(X(s]b,t),s)w(X(s]b,t),s)Xt(s]b,t)dsdt 

Jtn  Jtn 

ft*1*1  (i  —  tn)  / 

=  ~  Jta  - (/(M)  w(b,t)  Xt(t;b,t) 

+  f(X(tn]b,t),tn)iu(X(tn]b,t),tl)Xt(r]b,t))  dt+  Ef&(w) 
r^1  (/  /'■) 


(2.17) 


Jtn 


+ 


X 


At(0)  (a:) 


-1/ (6,  i)  f(b ,  f)  uj(6,  i)  dt 
f(x,tn )  r«(rc,t”)  da:  +  Ef^3(w), 


where  in  the  last  equality  we  used  the  second  equation  in  (2.10)  for  the  first  integral  and  the 
change  of  variable  x  =  X (I r>]  b.  I)  for  the  second.  The  trapezoidal  error  term  introduced  is 
given  by 


Ef.sh  M  =  /  +1  / 

,/t»+i  Jv 


ft  (tn  -  s)(t  -  s)  d2 


ds2 


[Xf(s;  6,  f)  f(X(s ;  6,  f),  s)]  ro(6,  f)  dsdt.  (2.18) 


2.2.2  Treatment  of  the  Diffusion  Term 


Next  we  treat  the  second  term  on  the  left  hand  side  of  equation  (2.4),  which  we  break,  similar 
to  the  source  term  in  (2.11),  over  the  three  subdomains  fl,  (i  =  1,2,3)  discussed  above.  Due 
to  the  similarity  in  the  treatment  of  this  integral  over  Qi  and  fi2,  we  combine  the  two  and 
use  the  same  change  of  variable  y  =  X(s\  aq£”+1)  as  before,  and  obtain 

//  D(y,  s)  uy(y,  s)  wy(y,  s)  dyds 

J J UO2 

=  [  I  Xx(s]  x,t,l+1)  (Dux)(X(s]  x,t,1+1),  s)  wx(X(s;x,tn+1),s)  dsdx 

Ja  ./•"  1..  i 

rb  /\f  / 

=  I  —  ([D^)(x,  r+1)  wx(x,  tn+1)  +  Xx(t*(x);  X,  tn+1)  x 

(Dux)(X(t*  (x)i  x,  I.'"  (x))  wx(X(t*(x)]  x,  i,,+1),  t*(x))j  dx 

f  A t(I)(x)  (Dux)(x,tn+1)  trx{x.  ')  dx  +  Ed^12(w) 

Ja 

A t{r)(x)  (Dux)(x,tn+1)  wx(x,tn+1)  dx  +  I  —  (Dux.)(x,tv+1)  wx(x,tn+1)  dx 

pb*  /\f 

+  /  —  {l)tix){.r.i")  /c,.;(.r.  O  dx+ 

./a  a 

(2.  i9) 

where  in  the  second  equality  we  used  a  backward  Euler  approximation  for  the  integral  over 
(a,  a)  and  a  trapezoidal  approximation  over  (a,  b )  with  error  given  below,  and  in  the  third 


+ 
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equality  we  made  the  change  w,,(X(s;  x,tn+1),  s )  =  wx(X(s\  x,tn+1 ),  s)x  Xx(s;  x,tn+1).  The 
error  term  in  (2.19)  is  given  by 


En,ntAw)  ~ 

ff“+1  rtn+1  d 


wx(x,  tn+1 )  dx 


+ 


I"  f  [  -(Dur)(X(6]x,tn+1)A)d6di 

Ja  Jt*(x)  Js  ad 

b  rt"+1  (fn  _  A(fT>  + 1  _  o')  d2  r  1 

- - - [(Dux)(X (s;  x,tn+1),  s)\  tr,.(x.t"~])  dsdx. 


(2.20) 


// 


Using  a  trapezoidal  approximation  for  this  diffusive  flux  term  over  (a,  a)  when  the  inflow 
boundary  condition  is  Diriclilet  or  Robin  condition  introduces  unknowns  on  the  inflow  bound¬ 
ary  and  severely  complicates  the  scheme.  Therefore  we  chose  to  use  only  backward  Euler 
approximation  on  that  interval.  Neumann  inflow  condition  on  the  other  hand  allows  us  to 
use  the  more  accurate  trapezoidal  approximation  for  the  diffusive  integral  uniformly  over 
the  whole  spatial  domain.  This  treatment  does  not  lead  to  any  non-symmetric  terms  in  the 
formulation.  Details  of  this  treatment  and  boundary  implementation  are  discussed  in  the 
following  section. 

Finally  we  treat  remaining  diffusion  integral  over  Q3.  With  the  change  of  variable  y  = 
X(s]b,t )  this  integral  becomes 


/ /  D(y,  s)  uy(y ,  s)  wy(y ,  s)  dyds 

=  ~  [  [  Xt(s-,b,t)  (Dux)(X(s;b,t),s)  wg(X(s;b,t),s)  dsdt 

Jtn  Jtn 


/•in+1  (l  _ 

Jta  — wt(bi b) 


(2.21) 


+Xt(f;b,  t )  ( Dux)(X(tn-b ,  t),  tn )  wx(X(tn ;  b,  t ),  tn+) )  dt  +  (w) 

ftn  +  1  U  _  fA 

=  ~  - x — ~(Dux)(b,t)  wt(b,t)  dt 

Jtn  Z 

rb  At(°Ux) 

+  - - - (Dux)(x,tn)  wx{^tn+)  dx  +  EDM3(w), 

where  in  the  second  equality,  we  have  used  the  relation  wx (b,t)Xt(t]b,t)  =  wt(b,t)  for  the 
first  integral.  The  error  term  U/j.o,,  is  given  by, 


f)(/  -  s)  dt2 
~2  d? 


( Dux)(X(s-,b,t ),  s)wt(X(s]  b,t),  s ) 


dsdt. 


(2.22) 
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2.2.3  Assembly  of  the  Reference  Equation 

Substituting  the  integrals  expanded  in  equations  (2.13),  (2.15),  (2.17),  (2.19)  and  (2.21)  back 
into  the  variational  equation  (2.4)  yields  the  following  reference  equation 


+ 

+ 


f  u(xyin+1)  w(x,tn+1)  dx  +  f  At{r)  (x)  (. Dux)(x,tn+1 )  wx(x,  tn+1)  dx 

J  a  J  a 

p  b  A  -t 

J  —  (t)iix)(xJ" '  ')  trjx.f  ')  dx 
pb  ( x)  ft n~^1  (f  —  t71') 

/  - t  -  ~  (l)ti,.)(.r.  I")  wx(x,t\)  dx  -  /  - - — -(I)  u,)(b.  I)  ir,  (I).  t. )  dt 

Ja  Z  Jtn  Z 

+  /  (Vu  —  Dux)(b,t)  w(b,t)  dt  —  /  (Vu  —  Dux)(a,t)  w(a,t)  dt 
Jtn  -ltn 

=  (  u(x,  tn)  w(x,  t")  dx  +  f  —  ^  ^ /(a,  t11+1)  w(x,  f1+1)  dx 

Ja  J a  2 

(' tn  +  1  ~  t) 


+ 


+ 


6  A/'"1  (x)  >‘tn+1  (+n+1 

— [-tf(x.r)  w[xj:)dx  + 


Jt 


V(a.  t)  f(a.  t)  ir  (a.  f.)  dt 


((-(") 


V (b,t)  f  (b,t)  w(b,t)  dt  —  E(w)  +  R„(w). 


(2.23) 

The  last  two  terms  in  equation  (2.23)  are  E(w)  =  YZl=\{E +  Eo.n ,)  which  represents  the 
collective  approximation  error  and  R„(w )  =  +/a6  u  (wt  +  Vwx)  dxdt  represents  the  adjoint 

term  which  in  our  case  vanishes  by  equation  (2.6). 


3  Numerical  Approximation 


Since  we  do  not  impose  a  particular  form  for  the  velocity  field  V(x,t),  explicitly  solving  the 
ordinary  differential  equation  (2.7)  which  defines  the  characteristics  is  not  possible  in  general 
and  introduces  additional  difficulty.  Therefore  in  our  numerical  approximation  of  the  solution 
of  the  advection  diffusion  equation  (2.1),  we  consider  an  approximation  of  the  characteristics 
X(8]x,t )  which  is  based  on  a  second  order  Runge-Kutta  approximation  known  as  Heun’s 
method.  In  particular,  we  define  the  approximate  characteristic  curve  emanating  from  a 
point  (x,t),  with  t  6  [tn,tn+1],  by 

y(0;  x,t)  =  x-  ZxZl  [V(x,t)  +  V(x  -  (t-  8)V(x,t),  0)] .  (3.1) 

Furthermore,  if  the  Courant  number  Or  =  max(V’) A// Ax  is  at  least  2,  we  partition  the 
outflow  boundary  {x  =  6,  f  G  [i”  ,Tl+1]}  as  follows: 


tt  = 


tn+1  -  (i  -  I)  Ax/  max(V),  i  =  /,...,/  +  IC  -  1, 

i  =  I  +  I C. 


(3.2) 


where  IC  in  this  case  is  the  integer  part  of  Cr.  When  no  subdivision  of  the  outflow  boundary 
is  introduced  (Cr  <  2),  we  let  IC  =  1  so  that  both  cases  are  treated  uniformly.  This  partition 
introduces  some  unknowns  on  the  outflow  boundary,  however  it  has  the  advantage  of  insuring 
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that  the  schemes  we  develop  are  suitable  even  for  large  Courant  numbers.  We  define  the 
test  functions  at  time  tn+1  as  hat  functions  given  by 


wt(x,t"+1)  = 

for  i  =  0, —  1.  At  the  outflow  boundary  the  test  functions  are  given  by, 

t  e  [t,_ i,t,:], 

te[ti,ti+ 1],  (3-4) 

otherwise, 

for  i  =  I  +  1, . . . ,/  +  IC.  Here  Ax;  =  X;  —  X;,_i  and  At;  =  t1_1  —  t,.  For  the  interior 
of  the  domain  Q,  we  extend  these  test  functions  defined  by  (3.3)  and  (3.4)  to  be  constant 
along  the  approximate  characteristics.  The  test  function  Wj  is  a  combination  of  both,  that 
is,  Wj(x,V1+1)  is  defined  by  (3.3)  for  the  interval  [xj_1,xr\,  and  wj(b,t )  is  defined  by  (3.4) 
on  the  interval  [f/+i,  tn+1]. 

The  numerical  schemes  are  based  on  approximating  the  exact  solution  u,  which  satisfies 
the  reference  equation  (2.23),  by  a  piecewise  linear  function  U  which  satisfies  a  similar 
equation  with  two  differences:  (i)  the  new  equation  will  be  developed  using  the  change  of 
variables  resulting  from  the  approximate  characteristic  tracking  given  by  (3.1),  and  (ii)  the 
error  and  adjoint  terms,  similar  to  E(w)  and  Rn  in  (2.23),  are  not  included.  Since  the 
terms  neglected  contribute  global  errors  which  are  of  order  0((Ax)2  +  (At)2),  the  resulting 
schemes  will  be  of  desired  accuracy.  In  the  following  two  subsections  we  develop  two  schemes 
to  solve  the  advection-diffusion  equation  (2.1)  based  on  backward  and  forward  characteristic 
tracking,  respectively.  Here  we  emphasize  that  the  test  functions  defined  above  are  constant 
along  the  approximate  characteristics. 

3.1  Backtracking  Runge-Kutta  characteristic  scheme 

The  first  numerical  scheme  (BRKC)  is  based  on  a  backward  tracking  of  characteristics.  We 
first  let  the  domains  ft,  be  defined  in  a  similar  manner  as  before,  but  use  the  approximate 
characteristics  Y(6\  x,ttn+1)  in  place  of  the  true  characteristics  (see  Figure  l.a).  In  this  case, 
our  characteristics  originate  at  points  ( x,tn+1 )  and  are  given  by  Y (6\x,V1+1)  in  fR  and 
0,2-  In  3 ,  the  characteristics  originate  at  points  (b,t)  on  the  outflow  boundary  and  are 
given  by  Y(6\b,t).  The  notation  introduced  in  Section  2  for  exact  tracking  is  also  used  in 
our  numerical  scheme  formulation,  however  we  modify  it  for  this  approximate  characteristic 
tracking.  Accordingly,  we  define  x  satisfying  x  =  Y (f";  x,  tn+1)  and  x*  =  Y(tn\  x,  tn+1)  as 
the  head  and  the  foot  respectively,  of  the  Runge-Kutta  approximate  characteristics.  The 
foot  of  a  characteristic  with  head  on  the  outflow  boundary  is  denoted  by  b*(t)  =  Y(tn ;  b,t) 
for  t  G  [tn,tn+1].  We  also  define  t*(x)  (given  by  the  relation  a  =  Y(t*(x );  x,t11+1))  as  the  exit 


V';(b.t.)  = 


R—  i  ^ 
At,:  ’ 
t  1 


At; 


*+l 


0, 


(  x  -  a,_i 
AX:t 

X'i-\- 1  X 


Axj 


*+i 


0, 


x  E  [aq_i,x,], 
x  G  [ x , , 
otherwise. 


(3.3) 
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time  of  the  characteristic  Y{6\  x.  t,1+1)  at  the  inflow  boundary.  On  the  other  hand  t ( x )  (given 
by  x  =  Y(tn]b,t(x)))  is  such  that  the  point  ( b,t(x ))  on  the  outflow  boundary  backtracks  to 
the  point  (. x,tn ).  The  two  time  increments  A t{T\x)  and  A t^(x)  are  as  defined  before  in 
(2.12). 

The  reference  equation  using  this  approximate  characteristic  tracking  is  derived  in  a 
similar  manner  as  was  done  for  equation  (2.23).  Since,  however,  the  solution  of  u  at  time 
level  t“+1  is  sensitive  to  errors  arising  from  the  evaluation  of  the  terms  of  (2.23)  which 
involve  integrals  of  the  trial  function  U(-,tn )  at  the  previous  time  level,  these  terms  must 
be  treated  with  care.  That  these  integrals  are  difficult  to  evaluate  (especially  in  higher 
dimensions)  is  due  to  the  fact  that  the  test  functions  are  defined  by  extending  back  along  the 
approximate  characteristics  their  values  from  time  tn+1  and  may  be  substantially  distorted 
by  this  process  [3,  45,  53].  To  indicate  the  possible  complications  that  may  arise,  we  note 
that  quite  complicated  geometries  may  occur  even  in  the  case  of  a  constant  velocity  field. 
For  example,  when  the  support  of  a  test  function  in  three  dimensions  is  traced  back  in  time 
and  intersects  a  boundary  of  the  spatial  domain,  a  four  dimensional  space-time  region  with 
a  complicated  geometry  results  as  the  support  of  the  space-time  test  function. 

The  backtracking  scheme  uses  the  approximate  characteristic  Y  to  change  variables  in 
each  of  the  integrals  where  the  test  functions  are  evaluated  at  time  level  f" .  In  this  way, 
the  integrals  are  rewritten  as  integrals  at  time  tn+1  (or,  as  the  case  may  be,  at  the  outflow 
boundary)  of  test  functions  (i.e.  hat  functions  on  the  grid)  at  the  current  time  integrated 
against  trial  functions  evaluated  at  the  foot  of  the  characteristics.  The  term  ‘backtracking’ 
comes  from  the  use  of  the  backward  characteristics  to  determine  the  trial  function  values  at 
the  previous  time.  Hence  after  some  rearrangement  and  combining  of  terms,  the  reference 
equation  of  the  piecewise  linear  trial  functions  U  using  Runge-Kutta  backward  tracking  is 
given  by 

f  U(x,  tn+1)  iu(x,  tn+1)  dx  4-  f  A/'  r'(.r)  ( t)l '..)(x.  t."  '  ')  tr  ..(x.  t"  '  ')  dx 

J  a  J  a 

pb  Af  , 

+  J  —  ((T>£4)(*0"+1)  +  (DV.:..)(x\t"))  wx(x,tn+1)  dx 

+  /  ( VU  -  DU,){b ,  t)  w(b ,  t)dt-  /  (VU  -  DUx)(a ,  t)  w(a,  t)  dt 

Jtn  Jtn 

-tn+i  (j-  _  +n\ 

-  J  {-^-L  (( DUx)(b,t )  +  (DUx*)(b*(t),tn)]  wt(b,t)  dt 

=  [b  Y^t^t**1)  U(x\tn)  w(x;f+1)  dx  -  f  Yit'>\b.t)U(h  (  t),t")  trih.l.)  dt. 

Ja  Jtn 

+  L  At  (/(M"+1)  +  Y.(f  (x):r.t"'  ')  f  (Y  ■  t:  (x)))  w(x,tn+1)dx 

«(»  +  !  ff  fn  \ 

+  l  (V(b,t)  J(b,t)  -  Y,(t";b,t)  J(b-[t),t’))w(b,t)  dt, 

'  ]  "  t  .  <3'5> 
where  in  the  third  integral  on  the  right  hand  side,  x*  is  a  on  (a,  a).  This  general  reference 

equation  holds  for  all  nodes  Xj,  {i  =  0, . . . ,  I)  and  i;,  (*  =  /,...,/  +  IC),  and  simplifies  in 

specific  cases.  The  stiffness  matrix  is  updated  by  evaluating  terms  involving  U  at  time  level 
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tn+1  and  at  the  outflow  boundary  in  equation  (3.5).  Likewise,  by  evaluating  the  remaining 
terms,  we  can  update  the  right  hand  side  vector  of  the  generated  system. 

As  we  mentioned  earlier,  one  difficulty  arises  in  accurately  computing  the  terms  evalu¬ 
ated  at  feet  of  characteristics  (i.e.  'ic*  or  since  several  grid  points  from  the  previous 

time  level  may  occur  as  backtracked  points  within  a  given  interval  at  the  current  time.  To 
illustrate  how  to  overcome  this  difficulty,  we  will  focus  on  the  first  term  on  the  right  hand 
side  of  equation  (3.5)  applied  with  test  function  w,.  In  this  case,  when  x  varies  over  the 
interval  [$&_!,  x-]  (i.e.,  the  left  half  of  the  support  of  m,(-,  fn+1)),  x *  will  vary  over  the  interval 
[•£;.  if  /'•]•  Therefore,  due  to  the  distortion  by  the  characteristic  tracking,  U(x*,tn ) 

is  not  a  piecewise  linear  function,  but  rather  a  piecewise  smooth  function  of  a:  £  [rc.,_i,  a|], 
This  possible  loss  of  smoothness,  however,  affects  the  accuracy  of  high  order  quadrature 
methods.  Therefore,  we  simply  subdivide  the  interval  so  that  U(x*,tn )  is  smooth  on  each 
subinterval  and  apply  numerical  integration  to  each  subinterval  separately.  To  illustrate  the 
idea,  we  suppose  x*_x  is  in  [.r .  , .  .r . )  and  x J  belongs  to  [sJ+i,  rcJ+2),  then  we  would  split  the 
integral  into  three  parts  along  the  intervals  (x.  | . .? ; ) .  [.? . .  J: . .  | ) .  and  [x,j+ i,s,).  Newton  iter¬ 
ation  is  applied  to  determine  the  points  at  the  current  time  level  that  track  back  to  the  grid 
points  at  the  previous  time  level.  (Recall  that  approximate  Runge-Kutta  forward  tracking 
and  back  tracking  may  not  be  inverse  operations  for  variable  velocity  fields).  Once  these 
intervals  are  determined,  we  may  numerically  integrate  by  determining  quadrature  points  in 
each  subinterval  and  then  backtracking  them  to  perform  the  quadrature  with  values  of  the 
trial  function.  The  amount  of  piecewise  smoothness  will  be  determined  by  the  corresponding 
smoothness  of  the  velocity  field.  Additional  information  on  this  general  problem  and  how  to 
treat  multivariate  problems  may  be  found  in  [28,  38,  52], 

Another  difficulty  in  incorporating  the  known  values  of  the  trial  function  in  equation  (3.5) 
arises  from  evaluating  the  expression  U^(x*,tv)  when  x*  happens  to  be  one  of  the  nodal 
points  {.x,}(_0  at  which  Ux*  may  have  a  jump  discontinuity.  Since  our  characteristic  tracking 
is  only  approximate,  we  assume  that  all  points  within  a  small  tolerance  of  x,  belong  to  the 
interval  (x,nx,i+1).  This  introduces  an  error  which  is  controlled  to  within  the  desired  order 
provided  the  tolerance  is  of  the  same  order. 

The  scheme  for  nodes  near  the  inflow  boundary  (i.e.,  x,  (i  =  0,  ...IC  +  1))  differs  since 
the  characteristics  traced  back  strike  that  boundary  and  lead  to  boundary  terms  appearing 
in  the  scheme.  In  the  case  of  Robin  flux  boundary  condition,  we  can  easily  change  the  fifth 
integral  on  the  left  side  of  equation  (3.5)  as  follows 

/  (VU  —  DUx)(a,t)  w(a,t)  dt  =  /  g3(t)  w(a,t)  dt.  (3.6) 

,/f»  Jt« 

Dirichlet  inflow  boundary  condition  introduces  extra  difficulty,  since  then  the  unknown  dif¬ 
fusive  flux  term  in  the  fifth  integral  on  the  left  side  of  equation  (3.5)  appears  in  the  scheme 
equations  for  nodes  x,  <  a.  An  alternative  treatment  for  the  diffusion  term  discussed  in 
[44,  52]  avoids  this  difficulty.  Instead  of  integrating  the  term  //„  —  (Dux)x  w  dxdt  by  parts 

as  we  did  in  the  derivation  of  the  variational  formulation  (2.4)  and  then  applying  the  integral 
approximation  (trapezoidal,  Euler),  we  can  reverse  the  order,  whence  we  get  the  term 

f  -  ^  ^  ( 1)1  V)(-r-  '  1 )  ir(xJ" '  1 )  dx  (3.7) 

J  a  dX 
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instead  of  the  diffusive  flux  in  the  fifth  term  on  the  left  side  of  (3.5).  The  fifth  term  on  the 
left  side  of  (3.5)  then  simplifies  to 

ftn+1 

/  V(a,t)  gi(t)  w(a,t)  dt.  (3.8) 

Jtn 

Since  the  value  of  the  trial  function  U  is  known  at  node  xq  from  the  prescribed  Diriclilet 
boundary  condition,  we  don’t  need  to  formulate  equations  there.  Hence  our  scheme  will  be 
stipulated  only  for  nodes  x,  (i  =  1 ,...,/).  However,  we  do  modify  w-\  =  w0  +  w1  so  that  the 
test  functions  sum  to  one  in  order  to  maintain  mass  conservation. 


Neumann  inflow  boundary  condition  generates  a  more  natural  scheme  in  the  sense  that 
a  trapezoidal  approximation,  instead  of  the  backward  Euler  approximation  described  above, 
can  be  implemented  for  the  diffusive  integral  over  fR.  The  advantage  of  this  treatment  is  that 
it  symmetrizes  the  inflow  terms  thus  maintaining  the  symmetry  of  the  whole  formulation. 
Implementing  this  boundary  condition  leads  to  the  following  two  changes  to  equation  (3.5): 
(i)  a  factor  of  1/2  multiplies  the  second  integral  on  the  left  side,  and  (ii)  the  additional  term 
of 

-f  At  ^ g2(t* (x))  wx(x,tn+1)  dx  (3.9) 

J  a  Z 

is  added  on  the  left  side.  The  fifth  term  on  the  left  side  of  equation  (3.5)  is  replaced  by 


/  (VU)(a,  t)  w(a,  t)  dt  +  /  $2^)  w(a,t)  dt 

■Jtn  Jtn 


(3.10) 


Wang,  Ewing,  and  Russell  [52]  describe  in  detail  several  possible  ways  to  treat  the  first 
integral  in  equation  (3.10).  Their  findings  suggest  combining  this  term  with  the  first  term 
on  the  left  hand  side  of  the  reference  equation  (3.5)  and  replacing  these  with  the  term 


/  U (a,  r  •  ')  w(x ,  tn+1)dx. 

Ja 


(3.11) 


Although  this  term  is  not  exact  for  a  variable  velocity  filed  V,  the  error  introduced  is  within 
the  desired  order. 


Due  to  the  fact  that  we  generate  unknowns  at  the  outflow  boundary  at  nodes  i(,  boundary 
treatment  here  becomes  very  important.  First  we  note  that  l  '(I),  t")  is  known  from  the 
previous  time  step  solution,  hence,  we  don’t  impose  an  equation  at  tI+Ic  =  t“ .  Therefore, 
as  we  did  earlier  in  the  case  of  Diriclilet  inflow  conditions,  we  redefine  the  test  function 
wi+jc-_ l  :=  tt)j+jc- 1  +  wj+jci  to  maintain  mass  balance.  For  test  functions  w-n  (?'  =  /  + 
1, +  IC  —  1),  the  terms  on  the  inflow  boundary  and  terms  evaluated  at  time  tn+1  of 
equation  (3.5)  vanish.  Thus  for  nodes  t,  (i  =  /  +  1, . . . ,  /  +  IC  —  1),  the  reference  equation 
(3.5)  simplifies  to  the  following 


/  (vr  ini  ){b.t)  w:(b.t)  ,ii 

Jtn 

-  Jta  (( DUr)(b,t )  +  (. DUx.)(b*(t),tn ))  wu(b,t)  dt 

=  ~  [t+  Yt(?;b ,  t)  U(b*(t),tn)  Wi(b ,  t)  dt 

Jtn 

/  I  _  in\ 

y-^-L  ( V(b,t )  f(b,t)  -  Y,(tn-,b,t)f(b-(t),r))  w,(b,t)  dt. 


(3.12) 


14 


We  describe  below  liow  this  equation  changes  for  all  three  types  of  outflow  boundary  condi¬ 
tions.  The  equation  for  node  xj,  has  in  addition  to  the  terms  in  (3.12)  some  other  terms  that 
are  in  equation  (3.5)  which  don’t  vanish  for  this  node.  Neumann  outflow  boundary  condition 
can  be  incorporated  in  equation  (3.12),  simply  by  making  the  change  (1)1  ,  ){h,  t)  =  h2(t ) 
in  the  first  and  second  integrals  on  the  left  hand  side  of  this  equation.  Similarly  the  flux 
outflow  condition  can  be  incorporated  in  equation  (3.12)  by  applying  the  condition  to  the 
first  integral  on  the  left  hand  side,  and  using  the  relation  — ( DUx)(b,t )  =  /13 (t)  —  (VU)(b,t) 
for  the  second  integral  on  the  left  hand  side  of  equation  (3.12).  In  the  case  of  outflow  Dirich- 
let  condition,  the  solution  U(b,in+1 )  is  known  from  the  boundary  condition,  therefore  the 
equations  for  i  <  I  —  1  decouple  and  can  be  solved  for  U (aq,  tn+1),  ( i  =  0 ...I  —  1).  Unless  we 
desire  the  values  of  the  derivative  of  the  solution  at  the  outflow  boundary,  we  don’t  generate 
any  equations  there. 

3.2  Forward  Tracking  Scheme 

In  this  brief  section  we  indicate  the  development  of  a  scheme  (FRKC)  based  on  a  ‘forward 
tracking’  which  is  a  feasible  alternative  for  the  backward  tracking  scheme  of  the  last  subsec¬ 
tion  especially  for  multidimensional  problems.  With  the  test  functions  as  defined  earlier  by 
equations  (3.3)-(3.4),  we  use  the  reference  equation  (2.23)  to  derive  the  numerical  scheme 
for  Runge-Kutta  forward  tracking.  Here  again  we  wish  to  avoid  the  expensive  task  of  back¬ 
tracking  the  geometry  to  perform  the  integration  of  the  test  functions  at  the  previous  time 
level.  Instead  we  perform  this  integration  by  using  the  fixed  spatial  grid  of  the  trial  func¬ 
tion  at  t  =  t"  and  apply  numerical  quadrature.  The  values  of  the  test  function  required 
by  the  quadrature  are  obtained  by  forward  tracking  the  quadrature  points  and  evaluating 
Wi(x,tn+1 )  or  its  derivatives  at  time  tn+1 .  This  forward  scheme  is  only  used  to  evaluate  cer¬ 
tain  terms  to  incorporate  known  values  and  therefore  does  not  lead  to  distorted  grids  that 
comprise  a  major  drawback  of  explicit  numerical  schemes. 

In  this  scheme,  the  characteristics  originate  at  points  (x,  tn)  and  are  given  by  Y(6;  x,  tn )  in 
0,2  and  O3  or  they  originate  at  points  (a,  t )  and  are  given  by  Y(6\  a ,  t )  in  fR.  Here  again  we  use 
the  same  notation  established  in  Sections  2  and  3,  except  we  redefine  them  for  this  particular 
forward  tracked  approximate  characteristics.  Accordingly,  we  define  x  =  Y(tn+1]  x,t”)  and 
x *  satisfying  x  =  Y (tn+1  \  x* \  tn)  as  the  head  and  foot  of  the  approximate  characteristic. 
Moreover,  we  let  t*(x)  (given  by  x  =  Y(tn+1]a,t*(x)))  and  t(x )  (given  by  b  =  )'(I(x):  x.  I,")) 
denote  the  exit  time  of  the  approximate  characteristics  at  the  inflow  and  outflow  boundaries, 
respectively.  The  time  increments  At11*  and  A R0)  are  then  given  by  equation  (2.12),  where 
t* |b)  and  t(x)  extend  to  t"  and  tn+1,  respectively,  if  they  are  defined  by  characteristics  in 
02.  The  reference  equation  for  the  forward  scheme,  derived  in  a  similar  manner  as  equation 
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+ 

+ 

+ 


(2.23),  is  given  by 

J  U(x,tn+1 )  w(x,tn+1 )  dx  +  J  A t(I]  (x)  (DUx)(x,tn+1)  wx(x,tn+1)  dx 

fb  Af 

I  —{DU,.){x.tn  ■  ')  trjx.t" '  ')  dx 

/  - r— -  ( DUx)(x ,  r)  ?«! (x.  t" )  dx  -  /  - - — -( DUx){b ,  t)  ^(6,  t)  dt 

Ja  2  Jtn  2 

ptn  + 1  /»£w  +  ! 

/  (Kt/  -  Dt4)(fe,  t)  w(6,  t)  dt  -  /  (Kt/  -  DUx)(a,  t)  w(a ,  t)  dt 

Jtn 

rt',+1  (tn+1  -  t ). 


j1  U(x,  tn )  to(a;,  f”  )  dx  +  j " 

,in+i  (  /  _  ) 


- V(a,t )  f(ci,t)  w(a,t )  dt 


y#(t;  fef|t),  t”)  / (6,  t)  ui(fe,t)  dt 


Jtn 

t  ——^-2^-f{x,  in+1)  w(x,tn+1)  dx  +  f'  Ati0)^  f(x,tn)  ,r(x.  C  )  dx, 
J a  2  .J  a  2  ' 


+ 


The  boundary  treatment  is  similar  to  that  of  the  backtracking  scheme  BRKC. 


(3.13) 


3.3  Experimental  Order  of  Convergence 

In  this  subsection  we  establish  numerically  the  order  of  convergence  of  the  two  schemes 
developed  for  equation  (2.1)  as  well  as  the  backward  Euler  ELLAM  [52],  referred  to  as  BE- 
ELLAM.  The  test  problem  involves  the  transport  of  a  Gaussian  distribution  (described  in 
Section  5.1)  initially  centered  at  zero.  The  spatial  domain  is  [0, 1.4]  with  a  temporal  domain 
[0, 1].  In  this  test  we  choose  the  small  diffusion  coefficient  of  D  =  10-4,  and  consider  the  two 
velocity  fields:  V  =  1  +  0.1s  and  the  more  rapidly  changing  V  =  1  —  0.5a:.  The  right  hand 
side  in  equation  (2.1)  is  generated  from  this  data  and  the  analytical  solution.  The  L2  and 
Li  norms  of  the  residual  error  of  the  solution  U  of  all  three  schemes  are  given  by 


max  II  U(x,tn) 

n=0,..,N  11  ' 


-  u(x,fl)\\LP(a,b]  <  Ca(Ax)a  +  Cl3(Aty 3, 


P=  1,2 


where  a,  (3  give  the  order  of  convergence  of  the  error  in  space  and  time  respectively,  and  Cr)s, 
C  j  are  positive  constants.  To  obtain  the  order  of  convergence  in  space  we  fix  At  =  1/500  to 
insure  time  truncation  errors  are  appropriately  small.  We  then  perform  runs  varying  Ax  and 
apply  a  linear  regression  on  the  Lp  (p  =  1,2)  norms  of  the  error  to  determine  the  parameter 
a.  Tables  1  and  3  present  the  results  of  these  runs  and  the  computed  value  of  the  parameter 
a  for  the  two  velocity  fields  described  above.  In  a  similar  manner  we  fix  Ax  =  1/700  and 
perform  runs  varying  At  to  estimate  the  value  of  the  parameter  (3.  Tables  2  and  4  give 
corresponding  results  for  (3  for  the  two  velocity  fields.  From  these  results  we  clearly  see  that 
the  two  schemes  developed  are  corroborated  to  be  second  order  in  space  and  time.  We  also 
see  that  BE-ELLAM,  as  was  expected,  is  second  order  in  space  but  only  first  order  in  time. 
The  orders  are  more  apparent  in  the  rapidly  changing  velocity  field  V  =  1  —  0.5.x  (Tables  3 
and  4).  Tables  1  and  3  show  that  when  the  time  step  At  is  very  small,  the  BRKC,  FRKC, 
and  BE-ELLAM  schemes  generate  solutions  with  comparable  errors.  This  is  expected  since 
all  the  schemes  are  second-order  accurate  in  space  and  the  temporal  errors  are  negligible.  In 
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practice  one  wishes  to  use  largest  possible  time  step  in  numerical  simulations  to  enhance  the 
efficiency  without  sacrificing  accuracy.  Therefore,  Tables  2  and  4  bear  more  computational 
aspect  since  very  large  time  steps  are  normally  desired.  One  sees  that  the  BRKC  and 
FRKC  schemes  further  reduce  the  temporal  errors  in  BE-ELLAM,  which  are  themselves 
significantly  smaller  than  most  Eulerian  methods.  This  justifies  the  appropriateness  of  these 
schemes  when  large  time  steps  are  to  be  taken.  The  constant  C p  is  much  smaller  than  Ca  for 
all  three  schemes  corroborating  that  time  truncation  errors  are  smaller  that  spatial  errors,  a 
strong  advantage  of  characteristic  methods  in  general. 


4  Description  of  Some  Other  Methods 

The  necessity  to  numerically  solve  advection-dominated  advection-diffusion  equation  with 
high  accuracy  has  lead  to  the  development  of  many  specialized  methods.  In  this  section 
we  briefly  describe  some  well  perceived  methods  which  are  widely  used  in  practice.  In  the 
next  section  we  carry  out  experiments  to  compare  the  performance  of  the  BRKC  and  FRKC 
schemes  with  these  methods.  In  this  description,  we  impose  all  the  assumptions  on  the 
velocity  field  V(x,t)  and  the  diffusion  coefficient  D(x,t )  which  were  described  in  Section  2 
and  describe  the  methods  with  Diriclilet  boundary  conditions. 


4.1  Galerkin  and  Petrov-Galerkin  finite  element  methods 


The  linear  Galerkin  (GAL),  quadratic  Petrov-Galerkin  (QPG),  and  cubic  Petrov-Galerkin 
(CPG)  FEM  methods  with  a  Crank-Nicliolson  time  discretization  and  Dirichlet  inflow  and 
outflow  boundary  conditions  for  equation  (2.1)  are  described  by  the  following  formulation 


pb  pb  /\f 

/  l  '(x.  I"  '  1  )(C;(x)<lx  /  -  (VU  —  1)1  ’....)  (X.  )  tl';,.(x)<lx 

J  a  J a  2 

pb  pb  A  -j 

=  /  l  '  (x.  t")(U:(x)dx  /  —  (VU  —  1)1  ’....)  (./’.  tr‘)tC;,.(x)<lx  (4-1) 

J  a  J a  2 

pb  A  f  ,  N 

+  J  ~2~  (/(M”+1)  +  /(#**"))  wt(x)dx 


where  the  trial  function  U ( x ,  tn+1 )  is  piecewise  linear  on  the  intervals  defined  by  the  partition. 
The  three  methods  differ  in  their  choices  of  the  test  functions  ui,{x)  (z  =  1,2,.../  —  1).  In 
the  GAL  method,  the  test  functions  are  chosen  from  the  space  of  the  trial  functions,  thus, 
w,(x)  are  the  hat  functions  defined  by  equation  (3.3)  where  w,  (x)  replaces  Wj(x,  t^1).  The 
Petrov-Galerkin  methods  were  designed  to  improve  the  GAL  method  by  introducing  some 
upwinding  in  the  test  space.  The  QPG  methods  use  test  functions,  which  are  the  sum  of  the 
standard  hat  functions  and  quadratic  asymmetric  perturbation  terms,  defined  by 


(  ( X  -  Xj_P)  +  3  ^  (x  -  xi_1)(xi  -  x) 


Wi(x )  =  < 


(te+i  -  x) 


Ax 


—  3  a 


H-i 


*+i  " 


( X  -  Xj)(x,l+ 1  -  x) 
Ax?+1 


x  e  p,_!,  aqj, 
x  e  [x.nx,  +  1], 
otherwise, 


(4.2) 
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where  the  parameter  a,  =  cotli( 1 )  —  with  V. ]  and  D,  being  the  arithmetic  mean  of 
the  velocity  field  and  the  diffusion  coefficient  over  the  interval  (x,_i,  xi).  The  CPG  methods 
on  the  other  hand  use  test  functions  which  have  symmetric  cubic  perturbation  terms  added 
to  the  hat  functions  as  follows, 


f  ( X  -  fl-i)  +  5  (X  -  -  x)(x.i_1  +  Xi  -  2x) 


W;(x)  =  l 


Ax, 

(x,i+1  -  x) 


Ax 


—  57; 


*+r 


*+i 


A-h! 

(x  —  Xj)(xj+ 1  —  x)(xi  +  x,i- |_i  —  2x) 
Axf+1 


{  o 


x  e  |.r.  I,  /;  . 

X  \  X , .  X  ,y  |  ] . 

otherwise, 


(4.3) 


where  7 |  =  CV,-2,  with  CV,:  =  AA-  js  the  Courant  number  averaged  over  the  interval  (xi_i,x.;). 
Although  many  other  Petrov-Galerkin  formulations  are  also  known  to  solve  equation  (2.1) 
reasonably  well,  QPG  and  CPG  methods  are  among  the  most  popular  ones  in  practice. 


4.2  The  Streamline  Diffusion  FEM  methods 


The  Streamline  Diffusion  method  (SDM)  is  applied  to  the  non-conservative  form  of  equation 
(2.1)  given  by 

Cu  =  ut  +  V(x,  t)ux  +  y:.(x.t)a  -  (/I (.r;,.f =  f(x,t).  (4.4) 


Here  we  describe  the  linear  SDM  formulation  for  this  equation.  We  divide  the  domain  into 
the  time  slabs  [a, b]  x  (tn,tn+1],  and  successively  on  each  slab,  we  seek  a  continuous  and 
piecewise  linear  function  U(x,t )  (discontinuous  in  time  at  t"  and  tn+1 )  which  satisfies  the 
following  formulation, 

j-tn  +  1  rb  ptn  +  1  rh 

/  /  Wt  +  VUX  +  VXU]  [ttJ  +  S(wt  +  Vwx)\  dxdt  +  /  /  DUX  wxdxdt 

Jt "  ./«  '  '  Jtn  Ja 

-S  f  f  ’\dUx)x  (w,  +  Vwx)  dxdt  +  [b  UyWydx  (4.5) 

Jtn  Ja  Ja 

ptnJt^  pb  pb 

=  If  [#+  S(wt  +  Vwx)\  dxdt  +  /  U^w^dx. 

JtnJa  Ja 


The  test  function  w  is  piecewise  linear  in  both  space  and  time  (bilinear)  in  the  slab,  and 
is  discontinuous  at  time  tn  and  V1+1  and  is  zero  at  x  =  a  and  x  =  b.  In  equation  (4.5) 
tUy  =  limi_qi>q±  w(x ,t),  Uf  =  u„(x),  and  S  is  a  free  parameter,  described  below,  that  has 
significant  influence  on  the  accuracy  of  the  scheme.  There  are  many  relations  that  can  be 
used  for  the  parameter  5.  One  of  the  most  widely  used  relations  is  <5  =  C— ==j,  where  h 
is  the  mesh  diameter  and  C  is  a  constant  to  be  chosen  [35].  The  choice  of  C  determines 
the  amount  of  diffusion  applied  in  the  direction  of  the  characteristics  and  therefore  has  a 
great  effect  on  the  accuracy  of  the  scheme.  One  requires  that  C  is  large  enough  to  produce 
non-oscillatory  solution,  but  not  too  large  to  damp  the  solution.  This  choice  is,  in  general, 
problem  dependent  and  not  clear  in  practice.  Although  there  are  more  improved  versions 
of  SDM  method  with  shock  capturing  capacity  which  produce  better  approximations,  they 
usually  have  nonlinear  formulations  (even  though  they  model  linear  equations),  have  more 
free  parameters  similar  to  <5  (described  above)  that  need  to  be  chosen  carefully,  and  also 
have  higher  computational  cost.  The  formulation  (4.5)  is  the  one  we  chose  for  numerical 
experiments  described  in  the  next  section. 
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4.3  The  Continuous  and  Discontinuous  Galerkin  FEM  methods 


The  Continuous  and  Discontinuous  Galerkin  methods  (CGM  and  DGM)  [42,  43]  apply  to  the 
non-conservative  form  of  equation  (2.1)  given  by  equation  (4.4).  The  domain  Qq  =  ( a,b )  x 
(0,T)  is  divided  into  a  quasi-uniform  triangulation  with  side  length  h,  and  Diriehlet  data  is 
assumed  on  the  inflow  portion  of  the  boundary  denoted  by  Tpj  and  given  by  ~V(x,t)  •  n  <  0, 
where  V (x.  t )  =  ( V ( x ,  t ),  1)  gives  the  characteristic  direction  and  n  =  (ni,  n 2)  is  the  outward 
unit  normal  vector.  The  boundary  data  given  at  the  outflow  portion  of  the  boundary  are 
not  used  in  the  formulation.  In  the  DGM  method  one  seeks  a  discontinuous  approximation 
U(x,t ),  which  lies  in  the  space  P„(T)  of  polynomials  of  degree  at  most  n  on  each  triangle 
T,  and  satisfies  the  equation 


+  VUr  +  VrU)  w  dxdt  -  f  U+  wV  ■  n  dS  +  [  DU+w  m  dS 

Jr  afiT)  •'r*J)(T) 

[  fw  dxdt  —  [  U~  w\  ■  11  dS  +  [  DU~w  rq  dS,  for  all  w  £  PAT ) 

Jt  Jr a)(T)  Jv;,  ,  /  i 

(4.6) 


where  T^(T)  is  the  inflow  boundary  of  T  exclusive  of  any  sides  on  the  boundary  of  the 
domain,  l ' '  (>•.  /)  =  liirn_>0+  U((x,t)  ±  eV),  and  U~  is  the  solution  at  the  previous  element 
or  an  interpolation  of  the  prescribed  Diriehlet  data  for  sides  on  T(/p  CGM  is  formulated  in 
the  same  way  by  requiring  a  solution  U(x,t )  6  P„(T)  which  satisfies  equat  ion  (4.6),  however 
there  are  two  main  differences:  First,  the  trial  function  U(x,t )  is  required  to  be  continuous 
over  the  domain  Qg,  and  the  test  functions  are  in  P„_p( t)(T),  where  p(T)  is  the  number  of 
inflow  sides  that  T  has.  Secondly,  the  continuity  requirement  in  CGM  makes  the  second 
terms  on  both  sides  of  (4.6)  cancel  each  other.  In  both  of  these  methods,  one  iterates  over 
the  elements  solving  a  linear  system  of  order  equal  to  the  degree  of  freedom  on  each  element 
which  makes  these  schemes  quasi-explicit.  In  the  next  section,  we  consider  the  lowest-order 
CGM  (n  =  2),  and  DGM  (n  =  1). 


4.4  High  Resolution  Methods  (MUSCL  and  ENO) 

High  resolution  methods  from  computational  fluid  dynamics  are  known  to  be  good  for  purely 
hyperbolic  equations.  An  extension  of  these  methods  to  equation  (2.1)  is  based  on  time- 
splitting  of  the  equation,  as  described  below,  and  then  a  higli-order  Godunov  method  can  be 
used  to  solve  the  advective  part,  with  a  mixed  FEM  method  to  solve  the  diffusive  part  [16]. 
We  consider  two  such  schemes,  the  first  based  on  Monotone  Upstream-centered  Scheme  for 
Conservation  Laws  (MUSCL)  which  was  developed  by  van  Leer  [49],  and  a  second,  based 
on  a  generalization  of  the  first,  called  the  Essentially  Non-Oscillatory  scheme  (ENO)  which 
was  developed  by  Harten  et  al.  [27,  46].  Assuming  that  Un(x )  approximates  the  solution 
u(x,  tn )  of  equation  (2.1),  we  can  generate  an  approximation  of  u(x,  f"+1)  as  follows:  First  the 
MUSCL  or  ENO  scheme,  described  below,  can  be  applied  to  find  a  solution  of  the  advective 
equation 

ut  +  (Vu)x  =  0,  for  (x,t)  £  [a,b\  x  (tv  ,tn+1], 

u(x,tn )  =  Uri(x),  for  x  £  [a,  b] , 
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which  we  denote  U,1+1(x).  Then  the  mixed  method  can  be  used  to  solve 


=  /, 

u*(x,tn)  =  Un+1(x), 


for  (x,t)  G  [a,  b]  x  (tn,tn+1], 
for  x  G  [a,  b] 


(4.8) 


whose  solution  is  the  approximation  U,1+1(x )  of  a(.t.  tn+1 ).  We  note  the  well-known  fact  that 
the  mixed  method  in  lowest-order  approximation  space  and  a  trapezoidal  rule  of  integration 
is  equivalent  to  the  block-centered  finite  difference  scheme  [21]. 

Now  we  describe  the  MUSCL  and  ENO  schemes.  Unlike  the  other  methods  discussed 
here  which  are  node  based,  MUSCL  and  ENO  are  cell-centered  based  methods;  i.e.  the 
solution  is  approximated  at  the  points  s,_  1/2  (*  =  1 where  £,-1/2  is  the  mid  point 
of  the  interval  To  simplify  the  presentation,  we  assume  that  the  partition  (2.2)  is 

uniform.  The  MUSCL  or  ENO  solution  of  equation  (4.7)  is  given  by 


At 


U"*%  =  u; ”_1/2  -  —  (v(n-t}aX)ui  1/2,L  -  v(xi_,x)uihL), 


(4.9) 


where  the  Courant  number  is  assumed  not  to  exceed  one  (a  stability  requirement  )  and  the 

'n 

i—l/2,L 


left  state  U"  /0  T  is  given  by 


Uh,2.L  =  < 


UI-1/2  +  —  _  V(Xi- 1/2,  *  —  1,  •  •  •  ,  /, 

i  9i(tn),  i  =  0. 


(4.10) 


The  two  methods  differ  in  the  choice  of  the  slope  SU™_1j2.  The  ENO  formulation  uses  the 
slope  given  by 


Wll/2  =  < 


^+Ui- 1/2: 


if  \A+U?_1/2\  <  |A _UV_ 


1/2  u 


A  U 


i- 1/25 


otherwise. 


(4.11) 


where  the  difference  operators  are  given  by 

I  2(tT1/2- 5i(C)) 


A  -Ullj2  =  < 


Ax 

TTn  —  77” 

L  i- 1/2  Ui- 3/2 

Ax 


i  =  1, 
i  =  2 


(4.12) 


and 


A +Ulm  =  < 


TTn  _  TTn 
Ui+ 1/2  C  /— 1/2 

As 

2(Mf)  -  c;i1/2) 

Ax 


i  =  l, 


(4.13) 


?;  =  /. 


The  MUSCL  formulation  on  the  other  hand  uses  the  following  definition  for  the  slope 

5Ullf2  =  min  {AlimU?_1/2,  |ACU"_1/2|}  •  sgn(AcX^1/2),  for  i  =  1, . . . ,  /,  (4.14) 


20 


where  sgn(s)  =  1  for  x  >  0,  sgn(x)  =  —1  for  x  <  1,  and  sgn(0)  =  0, 


'  a,min{|A+C/”1/2|,|A_C/t1/2|}, 

0, 


and 


AM, 


i- 1/2 


=  < 


Uj+ 1/2  -  4gi(tn) 
3Ax 

TTn  —  TJn 
L  i+ 1/2  L  i- 3/2 

2A®  ’ 


«  (A+EP  1/2)  •  (A_f/,”_1/2)  >  0, 


otherwise, 


(4.15) 


!  =  1. 

1  —  2, —  1,  (4.16) 


4h1(t’)-U’_3/2  . 

- : - — ,  *  =  I- 

3Ax 

The  parameter  at  in  equation  (4.15)  is  2  for  i  =  1,-  •  • ,  /  —  1  and  1  otherwise,  which  is  the 
upper  bound  that  allows  the  steeper  representation  of  sharp  fronts. 


5  Numerical  Experiments 

In  this  section  we  describe  numerical  experiments  which  we  use  to  compare  the  Runge-Kutta 
characteristic  methods  developed  in  this  paper  (using  both  forward-  and  back-tracking)  with 
several  generally  well-regarded  numerical  schemes,  including  various  Galerkin  and  Petrov- 
Galerkin  finite  element  methods,  high  resolution  methods  in  fluid  dynamics,  and  the  method 
of  Streamline  Diffusion.  We  apply  these  methods  to  two  standard  test  problems  (a  smooth 
Gaussian  distribution  and  a  step  function)  for  which  we  have  analytic  solutions  to  the 
advection-diffusion  equation.  In  addition,  each  of  these  functions  are  typically  used  to  test  for 
numerical  artifacts  of  proposed  schemes,  such  as  numerical  stability  and  diffusion,  spurious 
oscillations,  phase  errors,  and  Gibbs  type  effects  near  sharp  fronts. 

In  order  to  test  these  proposed  schemes  for  advecti on-dominated  transport,  we  consider 
the  model  equation  (2.1)  over  the  time  period  t  G  [0, 1.0]  with  a  variable  velocity  V(x,t)  = 
1+0.1  x,  and  a  relatively  small  diffusion  coefficient  of  D  =  10  4.  We  consider  both  test 
problems  for  our  initial/boundary  conditions  with  the  results  for  the  Gaussian  shown  in 
figures  beginning  with  label  I  and  the  step  function  with  label  II.  For  clarity  of  exposition, 
we  have  arranged  the  numerical  methods  into  5  groups  based  on  common  characteristics  of 
their  behavior  and  implementation.  These  groups  are  organized  according  to  the  following 
table: 


Group 

Methods 

1 

Runge-Kutta  characteristic  methods  (BRKC  and  FRKC) 

2 

Crank-Nicholson  FEM  (Galerkin,  Quadratic  and  Cubic  Petrov-Galerkin) 

3 

Streamline  Diffusion  (with  various  selections  of  the  control  parameter). 

4 

Continuous  and  Discontinuous  Galerkin 

5 

High  resolution  methods  in  fluid  dynamics  (MUSCL  and  ENO) 
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In  our  experiments,  we  have  systematically  varied  the  space  and  time  steps  to  examine 
the  performance  of  each  method.  For  each  grouping  we  have  chosen  to  display  3  plots 
which  provides  a  fair  illustration  of  the  accuracy  of  each  method,  their  potential  beneficial 
properties,  as  well  as  their  possible  undesirable  numerical  artifacts.  For  comparison  to  BRKC 
and  FRKC,  the  first  plot  in  each  figure,  labeled  (a),  presents  the  evolved  solution  for  all 
methods  in  the  group  using  a  common  space  mesh  (As  =  1/60  for  Problem  I,  Ax  =  1/100  for 
Problem  II)  and  a  time  step  as  close  as  possible  to  the  BRKC/FRKC  time  step  of  At  =  1/10, 
but  chosen  small  enough  to  ensure  stability.  The  third  plot  of  the  figure  (labeled  (c))  for  each 
grouping  shows  the  solutions  with  an  optimally  efficient  and  reasonable  choice  of  space  and 
time  steps  to  produce  a  qualitatively  comparable  solution  to  that  of  the  BRKC  and  FRKC 
schemes  in  Figures  1.1  and  II. 1,  respectively.  The  second  plot  in  each  Figure  (labeled  (b)) 
shows  an  intermediate  stage  for  each  grouping.  For  example,  Figure  1.2  consists  of  3  plots 
((a)-(c))  which  shows  the  solutions  for  model  problem  I  at  time  T  =  1.0  for  the  methods  in 
group  2  (spatial  finite  elements  and  a  Crank-Nicliolson  time  stepping)  with  (Ax,  At)  taken 
as  (1/60,1/69),  (1/60,1/200),  and  (1/100,1/300),  respectively.  In  this  example,  we  were 
required  by  the  CFL  constraint  to  begin  with  At  =  1/69,  since  the  maximum  velocity  is  1.14 
over  the  interval  [0,1.4]  and  larger  time  steps  result  in  unbounded  solutions. 

To  gauge  algorithm  efficiency,  we  also  compare  the  timings  for  each  method  to  achieve 
the  accuracy  depicted  in  plot  (c)  of  each  figure.  These  results  are  presented  in  Tables  5-9, 
using  the  groupings  of  algorithms  as  we  described  in  our  list  above.  We  realize,  of  course, 
that  some  code  optimization  may  be  possible  but  feel  that  these  timings  are  representative 
of  each  scheme’s  efficiency  on  these  model  problems. 

5.1  Model  Problem  I:  Gaussian 

Our  first  model  problem  is  a  Gaussian  distribution  over  the  spatial  domain  [0, 1.4].  We  use 
baseline  parameters  of  Ax  =  1/60,  At  =  1/10  for  each  numerical  scheme  in  our  testbed  and 
vary  these  parameters  until  we  obtain  results  with  errors  comparable  to  the  two  Runge-Kutta 
methods  (BRKC  and  FRKC).  With  the  exception  of  the  BRKC  and  FRKC  schemes,  each 
of  these  methods  has  implicit  Courant  restrictions  on  the  time  step  which  must  be  met  for 
numerical  solutions  to  be  bounded  and  may  not  be  chosen  as  large  as  that  permitted  by  the 
ELL  AM  schemes. 

Figure  1.1  shows  the  initial  condition  for  model  problem  I  (plotted  with  solid  line  at  the 
inflow  boundary)  along  with  the  evolved  solution  at  time  t  =  1.0.  Near  the  right  boundary,  we 
see  the  analytic  solution  (solid  line),  the  back-tracked  BRKC  solution  (marker  ’o’),  and  the 
forward-tracked  FRKC  solution  (dotted  line).  Both  methods  use  a  spatial  step  of  Ax  =  1/60 
with  a  relatively  large  time  step  of  At  =  1/10.  Each  of  the  two  solutions  give  very  accurate 
approximations  which  are  free  of  numerical  oscillations,  artificial  diffusion,  phase  error,  and 
adverse  boundary  effects.  The  timings  for  the  BRKC  and  FRKC  schemes  are  presented  in 
Table  5  and  provide  baseline  timings  for  all  experiments. 

Figure  1.2. a  contains  the  plots  of  the  analytic  solution  (solid  line),  the  Galerkin  ap¬ 
proximation  (labeled  with  symbol  ’+’),  quadratic  Petrov-Galerkin  (dotted  line)  and  cubic 
Petrov-Galerkin  (labeled  with  symbol  ’o’)  at  time  t  =  1.0  with  Ax  =  1/60  where  the  time¬ 
stepping  method  employed  is  Crank-Nicliolson.  Due  to  Courant  number  restrictions  on  these 
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methods,  an  initial  time  step  of  At  =  1/69  was  required.  This  plot  shows  that  there  are 
significant  trailing  oscillations  for  both  the  Galerkin  and  quadratic  Petrov-Galerkin  methods 
and  to  a  lesser  extent  for  the  cubic  Petrov-Galerkin  method.  All  methods  in  this  group 
however  have  significant  diffusion  and  a  mild  downstream  phase  error.  As  we  decrease  our 
mesh  sizes  to  try  to  match  the  performance  of  the  two  Runge-Kutta  characteristic  schemes, 
we  see  in  Figures  I.2.b-I.2.c  that  the  trailing  oscillations  and  numerical  diffusion  are  less 
pronounced  (for  all  but  the  cubic  Petrov-Galerkin  method),  but  still  persist  for  all  three 
methods  until  we  decrease  Ax  to  1/100  and  At  to  1/300.  In  this  case  the  CPU  requirement 
for  these  methods  is  two  orders  of  magnitude  larger  than  that  required  to  achieve  similar 
results  using  the  BRKC  or  FRKC  methods  (see  Table  6).  In  displaying  the  timings  for  these 
Crank-Nicholson  schemes,  we  have  presented,  in  the  interest  of  space,  only  the  average  of 
the  three  schemes  since  the  timings  for  each  scheme  was  within  one  second  of  the  other  two 
methods. 

Figure  1.3  presents  the  corresponding  results  for  the  Streamline  Diffusion  Method.  This 
method  requires  the  use  of  a  control  parameter  C,  and  we  present  in  Figure  1.3  the  plots  for 
three  values  of  this  parameter  in  each  of  the  subplots  (a)-(c).  In  Figure  1. 3. a  we  have  used  a 
time  step  of  At  =  1/120  which  is  the  largest  we  may  choose  and  have  bounded  solutions  due 
to  the  method’s  Courant  constraint.  As  Figure  1.3  demonstrates,  there  are  both  leading  and 
trailing  oscillations,  along  with  relatively  strong  numerical  diffusion  and  a  downstream  phase 
error.  These  persist  to  a  milder  degree  (see  Figure  1.3. c)  as  both  the  spatial  and  time  steps 
are  decreased.  In  Table  7  we  present  the  timings  for  this  method.  For  each  selection  of  C . 
the  CPU  time  is  presented  to  compute  the  solution  for  t  =  1.0.  The  timings  for  Figure  I.3.a-c 
are  at  least  68.7,  269.1,  and  402.1  seconds,  respectively,  which  indicate  a  severe  weakness 
of  this  method  in  certain  applications.  Another  drawback  of  this  method  is  that  it  is  not 
clear  in  general  how  to  choose  the  parameter  C,  which  indicates  that  an  iteration  of  this 
parameter  may  be  necessary  and  the  expense  of  this  method  will  increase  accordingly. 

Figure  1.4  presents  the  results  for  the  Continuous  (dotted  line)  and  Discontinuous  (labeled 
with  symbol  ’+’)  Galerkin  methods.  Again  due  to  Courant  restrictions,  we  were  able  to 
take  At  no  larger  than  Ax  =  1/60.  We  see  that  there  are  significant  leading  and  trailing 
oscillations  and  a  mild  phase  error  for  these  two  schemes  with  the  Discontinuous  Galerkin 
method  performing  somewhat  better  of  the  two.  However,  this  later  method  does  exhibit 
some  overshoot  near  the  maximum  of  the  Gaussian,  which  persists  in  plot  (b)  where  Ax  = 
1/120  and  At  =  1/180.  The  accuracy  of  the  two  Runge-Kutta  characteristic  schemes  are 
matched  in  plot  (c)  by  taking  Ax  =  1/180  and  At  =  1/180,  but  both  the  CGM  and  DGM 
require  a  CPU  expense  of  more  than  120  seconds  (see  Table  8)  as  compared  to  1.1  seconds 
for  the  FRKC  and  2.1  seconds  for  the  BRKC  (see  Table  5). 

The  final  schemes  which  we  wish  to  observe  are  the  two  high  resolution  methods  in  fluid 
dynamics  (MUSCL  and  ENO)  In  Figure  1.5,  we  have  plotted  the  results  of  the  simulation 
of  these  schemes  where  we  again  must  take  a  relatively  small  initial  time  step  of  At  =  1/69 
due  to  the  Courant  restriction  of  these  methods.  The  analytic  solution  is  plotted  as  the 
solid  line,  while  the  ENO  scheme  uses  a  dotted  line  and  the  MUSCL  scheme  is  plotted 
using  the  symbol  ’+’.  The  monotonicity  of  these  methods  is  quite  apparent,  but  there 
is  a  pronounced  trailing  non-negative  oscillation  as  well  as  an  overshoot  near  the  peak  of 
the  pulse  for  both  methods,  with  the  ENO  scheme  producing  the  most  overshoot.  These 
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effects  indicate  strong  mass  balance  errors  for  the  Godunov  schemes.  Parts  (b)  and  (c)  of 
this  figure,  show  that  these  artifacts  persist  until  we  reduce  the  mesh  sizes  to  have  values 
Ax  =  1/300  and  At  =  1/500.  Even  with  this  reduction  the  trailing  nonnegative  oscillation 
is  still  apparent  for  both  methods  (see  Figure  1.5. c).  To  achieve  this  level  of  approximation 
required  a  total  of  30.8  and  32.7  seconds  of  CPU  time  for  the  MUSCL  and  ENO  schemes, 
respectively,  which  is  a  significant  increase  over  the  1.1  sec  required  for  the  FRKC  scheme, 
or  for  the  BRKC  (2.1  seconds). 


5.2  Model  Problem  II:  Step  Function 

In  this  section  we  discuss  the  corresponding  experiments  carried  out  above  with  the  step 
function  replacing  the  Gaussian  distribution.  The  observations  made  from  Problem  I  are 
still  valid  for  all  of  the  numerical  schemes  tested,  but  in  this  case  there  are  several  additional 
features  we  would  like  to  point  out  which  are  less  pronounced  than  in  the  case  of  a  Gaussian 
distribution. 

For  this  model  problem  our  interval  is  [0,2]  and  we  begin  with  Ax  =  1/100.  Figure  II.  1 
shows  the  initial  condition  and  analytic  solution  (solid  line),  together  with  the  BRKC  and 
FRKC  solutions.  We  have  taken  At  =  1/10  for  these  methods. 

The  Galerkin  and  Petrov-Galerkin  methods  with  Crank-Nicholson  time  stepping  (see 
Figure  II. 2)  shows  strong  oscillations  (both  up-  and  down-stream)  near  each  steep  front. 
Figure  II. 3  shows  that  the  streamline  diffusion  method  exhibits  a  Gibbs’  effect  near  each 
of  the  jump  discontinuities  in  the  initial  function.  In  this  case,  it  appears  that  reducing 
the  value  of  C  continues  to  improve  the  error,  while  in  previous  experiments  [54],  we  have 
shown  that  the  error  may  have  local  minimum  as  a  function  of  C  for  this  method.  The 
Continuous  and  Discontinuous  Galerkin  methods  (see  Figures  II. 4  parts  b-c)  have  a  behavior 
very  similar  to  the  Streamline  Diffusion  method.  This  is  to  be  expected  since  these  methods 
are  closely  related.  However,  the  Streamline  Diffusion  method  has  approximately  double  the 
unknowns  which  must  be  determined  and  is  correspondingly  more  expensive  to  compute. 
The  high  resolution  schemes  (Figure  II. 5)  tend  to  worsen  slightly  as  the  grids  are  refined 
past  Ax  =  1/100  and  At  =  1/200.  Although  the  numerical  error  is  somewhat  high,  the 
qualitative  features  of  these  solutions  are  good. 
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Method 

At 

Ax 

1/2  Error 

L\  Error 

BRKC 

1/500 

1/40 

1.708301  x  10“2 

1.008468x  10-2 

1/500 

1/50 

9.762816  x  10“3 

5.012535x  10-3 

1/500 

1/60 

6.028822  x  10“3 

2.822963x  10-3 

1/500 

1/70 

3.989709  x  10“3 

1.738832x  10-3 

1/500 

1/80 

2.793234  x  10“3 

1.191840x  10-3 

a  =  2.6181 

a  =  3.0979 

=  270.4686 

Ca  =  919.3184 

FRKC 

1/500 

1/40 

3.280644  x  10“2 

1. 560330  xl0“2 

1/500 

1/50 

2.993149  x  10“2 

1.335517x  10-2 

1/500 

1/60 

2.290385  x  10“2 

1.005026x  10-2 

1/500 

1/70 

1.415099  x  10“2 

6.052967  x  10-3 

1/500 

1/80 

5.040757  x  10“3 

2.071540x  10-3 

a  =  2.4812 

a  =  2.6757 

=  418.3782 

Ca  =  405.2484 

BE-ELLAM 

1/500 

1/40 

1.727842  x  10“2 

1.019763x  10-2 

1/500 

1/50 

1.004383  x  10“2 

5.125955x  10-3 

1/500 

1/60 

6.370859  x  10“3 

2.974945  xl0“3 

1/500 

1/70 

4.375513  x  10“3 

1.910709x  10-3 

1/500 

1/80 

3.212402  x  10“3 

1.356822x  10-3 

a  =  2.4374 

a  =  2.9236 

=  138.478 

Ca  =  481.5463 

Table  1:  Order  of  convergence  in  space  with  V(x,t)  =  1  +  Q.\x 
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Method 

At 

Ax 

1/2  Error 

L\  Error 

BRKC 

1/6 

1/700 

5.653338  x  lO”3 

2.079667  x  10“3 

1/8 

1/700 

3.189920  x  10“3 

1.225680  x  10-3 

1/10 

1/700 

2.066911  x  10“3 

8.268313  x  10“4 

1/12 

1/700 

1.460440  x  10-3 

5.999894  x  10“4 

1/14 

1/700 

1.094716  x  10“3 

4.562014  x  10“4 

P  =  1.9385 

P  =  1.7867 

6%  =  0.1808 

C.,  =  0.0508 

FRKC 

1/6 

1/700 

5.920335  x  10“3 

2.503178  x  10-3 

1/8 

1/700 

3.310857  x  10“3 

1.400726  x  10-3 

1/10 

1/700 

2.117405  x  10-3 

8.985797  x  10“4 

1/12 

1/700 

1.467171  x  10“3 

6.138583  x  10“4 

1/14 

1/700 

1.075166  x  10-3 

4.520550  x  10“4 

(3  =  2.0123 

P  =  2.0222 

Cj  =  0.2177 

Cj  =  0.0939 

BE-ELLAM 

1/6 

1/700 

6.476122  x  10“2 

2.883357  x  10“2 

1/8 

1/700 

4.817777  x  10“2 

2.151369  x  10“2 

1/10 

1/700 

3.836627  x  10“2 

1.716823  x  10“2 

1/12 

1/700 

3.187941  x  10“2 

1.428786  x  10“2 

1/14 

1/700 

2.727120  x  10“2 

1.223664  x  10“2 

(3  =  1.0207 

3  =  1.0115 

Cj  =  0.4028 

Cj  =  0.1764 

Table  2:  Order  of  convergence  in  time  with  V(x,t)  =  1  +  0.1a; 
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Method 

At 

Ax 

1/2  Error 

L\  Error 

BRKC 

1/500 

1/60 

6.117323  x  10“3 

2.475200  xlO”3 

1/500 

1/70 

3.495166  x  10“3 

1.323890x  10-3 

1/500 

1/80 

2.429728  x  10“3 

9.213634x  10-4 

1/500 

1/90 

1.847878  x  10“3 

6.923624x  10-4 

1/500 

1/100 

1.408213  x  10“3 

5. 172978  xl0“4 

a  =  2.8265 

a  =  2.9903 

CQ  =  610.3031 

CQ  =  475.0919 

FRKC 

1/500 

1/60 

1.159784  x  10“2 

5.085293  xlO-3 

1/500 

1/70 

8.458101  x  10“3 

3.559665x  10-3 

1/500 

1/80 

5.902566  x  10“3 

2.528298  xlO”3 

1/500 

1/90 

3.967780  x  10“3 

1.713444x  10-3 

1/500 

1/100 

3.488353  x  10“3 

1.511827x10“  3 

a  =  2.4833 

a  =  2.4856 

CQ  =  308.5422 

=  134.2145 

BE-ELLAM 

1/500 

1/60 

6.673372  x  10“3 

2.733189x  10“3 

1/500 

1/70 

4.124203  x  10“3 

1.573640  xl0“3 

1/500 

1/80 

3.124454  x  10“3 

1.199037  x  10“' 3 

1/500 

1/90 

2.606328  x  10“3 

9.991517x  10“4 

1/500 

1/100 

2.240822  x  10“3 

8.839847  x  10“4 

a  =  2.1018 

a  =  2.1665 

Ca  =  33.5250 

=  17.3620 

Table  3:  Order  of  convergence  in  space  with  V(x,t)  =  1  —  0.5a: 
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Method 

At 

Ax 

1/2  Error 

Li  Error 

BRKC 

1/6 

1/700 

3.971642  x  10“2 

1.457259  x  10“2 

1/8 

1/700 

2.157900  x  10“2 

8.290132  x  lO”3 

1/10 

1/700 

1.387359  x  10“2 

5.448085  x  10“3 

1/12 

1/700 

9.796209  x  10“3 

3.888443  x  10“3 

1/14 

1/700 

7.347346  x  10“3 

2.936484  x  10“3 

(3  =  1.9897 

3  =  1.8896 

=  1.3770 

Cp  =  0.4261 

FRKC 

1/6 

1/700 

3.866812  x  10“2 

1.336697  x  10“2 

1/8 

1/700 

2.044585  x  10“2 

7.224535  x  10“3 

1/10 

1/700 

1.276627  x  10“2 

4.564537  x  10“3 

1/12 

1/700 

8.764417  x  10“3 

3.148085  x  lO”3 

1/14 

1/700 

6.401563  x  10“3 

2.313990  x  lO”3 

(3  =  2.1208 

3  =  2.0692 

Cp  =  1.7050 

Cp  =  0.5394 

BE-ELLAM 

1/6 

1/700 

1.267779  x  lO”1 

5.017525  x  10“2 

1/8 

1/700 

9.377586  x  10“2 

3.749495  x  10“2 

1/10 

1/700 

7.465792  x  10“2 

3.004816  x  10“2 

1/12 

1/700 

6.210892  x  10“2 

2.512159  x  10“2 

1/14 

1/700 

5.320968  x  10“2 

2.160284  x  10“2 

(3  =  1.0242 

3  =  0.9941 

Cp  =  0.7916 

Cp  =  0.2971 

Table  4:  Order  of  convergence  in  time  with  V ( x ,  t)  =  1  —  0.5a: 
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Method 

At 

Ax 

Z/2  Error 

L\  Error 

CPU 

Fig. 

BRKC 

1/5 

1/60 

8.370973  x  lO”03 

3.151613  x  10“03 

1.2 

- 

1/10 

1/60 

2.581547  x  10-03 

9.348955  x  10“04 

2.1 

1.1 

FRKC 

1/5 

1/60 

8.426960  x  10“03 

3.571262  x  10“03 

0.7 

- 

1/10 

1/60 

2.628836  x  lO”03 

1.039830  x  10“03 

1.1 

1.1 

Table  5:  Results  for  BRKC  &  FRKC.  CPU  is  in  seconds. 


Method 

At 

Ax 

Lo  Error 

L\  Error 

CPU 

Fig. 

GAL 

1/69 

1/60 

6.315931  xl0“2 

3.769911  xl0“2 

17.9 

I.2a 

1/200 

1/60 

1.384596  xl0“2 

6.736512  xlO-3 

51.8 

I.2b 

1/300 

1/60 

8.437177  x  10“ 3 

4.054622  xlO-3 

77.6 

- 

1/115 

1/100 

3.011717  x  10-2 

1.449988  xl0“2 

49.7 

- 

1/200 

1/100 

1. 159494  xlO”2 

5.076671  xlO”3 

86.5 

- 

1/300 

1/100 

5.493178  xlO”3 

2.319413  xlO”3 

129.2 

I.2c 

QPG 

1/69 

1/60 

5.042626  x  10-2 

2.597113xl0-2 

17.9 

I.2a 

1/200 

1/60 

2.384423  x  10-2 

1.121906xl0-2 

51.9 

I.2b 

1/300 

1/60 

2.596181  x  10“  2 

1.225615xl0-2 

77.7 

- 

1/115 

1/100 

2.564857  xl0“2 

1.212395xl0-2 

49.8 

- 

1/200 

1/100 

9.187296x  10-3 

4.064321  xlO”3 

86.5 

- 

1/300 

1/100 

6.861496  xl0“3 

2.966565  xlO-3 

129.3 

I.2c 

CPG 

1/69 

1/60 

1.502762  xl0“2 

6.299219  xlO-3 

18.0 

I.2a 

1/200 

1/60 

4.935717  x  10-3 

2.350595  xlO-3 

51.9 

I.2b 

1/300 

1/60 

4.589138  x  10-3 

2.242219  xlO”3 

77.8 

- 

1/115 

1/100 

5.494951  x  10-3 

2.272993  xlO”3 

50.0 

- 

1/200 

1/100 

1.897574x  10-3 

7.277279  xlO“4 

86.3 

- 

1/300 

1/100 

1.080232  xl0“3 

4.393256  xlO“4 

129.3 

I.2c 

Table  6:  Results  for  GAL,  QPG,  &  CPG.  CPU  time  is  in  seconds. 
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c 

At 

Ax 

L2  Error 

Li  Error 

CPU 

Fig. 

1.0 

1/60 

1/60 

1.123311  x  10-1 

6.024287  x  10-2 

68.7 

1. 3a 

1/120 

1/60 

8.824366  xlO-2 

4.589660  xlO”2 

134.3 

- 

1/180 

1/60 

7.944206  xl0“2 

4.087597  x  10-2 

197.9 

- 

1/60 

1/120 

8.795492  xlO-2 

4.613883xl0-2 

141.7 

- 

1/120 

1/120 

5.112340xl0-2 

2.561269xl0-2 

269.1 

1. 3b 

1/180 

1/120 

3.812530  xl0“2 

1.867242  xl0“2 

402.1 

1. 3c 

0.1 

1/60 

1/60 

5.241310xl0-2 

2.723927  x  10-2 

68.6 

1. 3a 

1/120 

1/60 

3.184598xl0“2 

1.581704x  10-2 

133.6 

- 

1/180 

1/60 

2.554619  xlO”2 

1.238932  xlO”2 

198.2 

- 

1/60 

1/120 

3.683633  xl0“2 

1.838774x  10-2 

140.4 

- 

1/120 

1/120 

1.514924x  10-2 

6.961759xl0-3 

270.4 

1. 3b 

1/180 

1/120 

9.512487  x  10“  3 

4.181123X10-3 

402.6 

1. 3c 

0.0001 

1/60 

1/60 

3.279756  xl0“2 

1.715544x  10-2 

69.1 

1. 3a 

1/120 

1/60 

1.623756  xl0“2 

8.350265  xl0“3 

140.4 

- 

1/180 

1/60 

1.156853  xl0“2 

5.913949  xl0“3 

201.2 

- 

1/60 

1/120 

2.310494x  10-2 

1.108312xl0-2 

149.7 

- 

1/120 

1/120 

7.797575  xl0“3 

3.349920  xlO”3 

284.2 

1. 3b 

1/180 

1/120 

4.380929  xl0“3 

1.731336xl0-3 

412.4 

1. 3c 

Table  7:  Results  for  SDM  with  C=l,  0.1,  &  0.0001,  CPU  is  in  seconds. 
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Method 

At 

As 

L-2  Error 

L\  Error 

CPU 

Fig. 

CGM 

1/60 

1/60 

4.574094  x  10“2 

2.216362x  10-2 

13.9 

1. 4a 

1/120 

1/60 

3.136981  x  10“2 

1.399978  xl0“2 

27.7 

- 

1/180 

1/60 

2.674044  x  10“2 

1.150175x  10-2 

41.5 

- 

1/60 

1/120 

2.501450  x  10“2 

1.191302x  10-2 

28.4 

- 

1/120 

1/120 

1.025814  x  10“2 

4.547582  x  10-3 

55.3 

- 

1/180 

1/120 

6.763004  x  10-3 

2.919515x  10-3 

85.1 

1. 4b 

1/60 

1/180 

1.956744  x  10“2 

9.230535  xl0“3 

41.5 

- 

1/120 

1/180 

6.110354  x  10-3 

2.687964x  10-3 

82.9 

- 

1/180 

1/180 

3.381005  x  10-3 

1.457474x  10-3 

124.3 

1. 4c 

CGM 

1/60 

1/60 

3.418172  x  10“2 

1.749322  xl0“2 

18.1 

1. 4a 

1/120 

1/60 

2.197196  x  10“2 

1.051318x  10-2 

36.2 

- 

1/180 

1/60 

1.810357  x  10“2 

8.461050x  10-3 

54.3 

- 

1/60 

1/120 

2.206570  x  10“2 

1.083687  x  10-2 

36.2 

- 

1/120 

1/120 

9.206862  x  10“3 

4.198669x  10-3 

72.3 

- 

1/180 

1/120 

6.763004  x  10“3 

2.727802  xl0“3 

85.1 

1. 4b 

1/60 

1/180 

1.855220  x  10“2 

8.954882  x  10-3 

54.3 

- 

1/120 

1/180 

6.062851  x  10“3 

2.722088  x  10-3 

108.5 

- 

1/180 

1/180 

3.502214  x  10“3 

1.536899  xl0“3 

162.7 

1. 4c 

Table  8:  Results  for  CGM  &  DGM,  CPU  is  in  seconds. 
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Method 

At 

Ax 

L2  Error 

L\  Error 

CPU 

Fig. 

MUSCL 

1/69 

1/60 

5.131587  x  10“2 

1.870076  x  10“2 

0.9 

1. 5a 

1/115 

1/100 

3.019053  x  10“2 

1.019578  x  10“2 

2.6 

1. 5b 

1/300 

1/100 

2.737727  x  10“2 

1.161478  x  10“2 

6.1 

- 

1/229 

1/200 

1.327108  x  10“2 

4.449515  x  10-3 

9.4 

- 

1/300 

1/200 

1.203944  x  10“2 

4.042067  x  10-3 

12.4 

- 

1/343 

1/300 

7.094933  x  10“3 

2.403457  x  10“3 

21.6 

- 

1/400 

1/300 

6.580204  x  10“3 

2.243087  x  10“3 

24.8 

- 

1/500 

1/300 

5.879445  x  10“3 

1.997493  x  10-3 

30.8 

1. 5c 

1/700 

1/400 

3.108750  x  10“3 

1.169307  x  10-3 

58.2 

- 

ENO 

1/69 

1/60 

5.096013  x  10“2 

1.966867  x  10“2 

0.9 

1. 5a 

1/115 

1/100 

3.042145  x  10“2 

9.933304  x  lO”3 

2.6 

1. 5b 

1/300 

1/100 

4.772205  x  10“2 

2.364377  x  10“2 

6.6 

- 

1/229 

1/200 

1.667658  x  10“2 

6.502723  x  10“3 

10.1 

- 

1/300 

1/200 

1.386139  x  10“2 

5.156693  x  10-3 

13.2 

- 

1/343 

1/300 

1.343450  x  10“2 

5.318024  x  10“3 

22.5 

- 

1/400 

1/300 

1.251113  x  10“2 

4.970423  x  10“3 

27.1 

- 

1/500 

1/300 

1.094915  x  10“2 

4.439829  x  10“3 

32.7 

1. 5c 

1/700 

1/400 

1.090687  x  10“2 

4.150278  x  10“3 

61.7 

- 

Table  9:  Results  for  MUSCL  &  ENO,  CPU  is  in  seconds. 
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anlaytic, "+"  GMM-MUSCL, GMM-ENO,dx=1/60,  dt=1/69 


(a)  Ax  —  60,  A t  —  6g 


anlaytic, "+"  GMM-MUSCL, GMM-ENO,  dx=1/100,  dt=1/1 15 


(b)  Ax  —  1J0,  At  — 


anlaytic, "+"  GMM-MUSCL, GMM-ENO,  dx=1/300,  dt=1/500 


( c )  A  =  — - —  A  f  =  — - — 

>  300’  500 


Fig  1.5  High  resolution  methods  MUSCL  and  ENO 
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1.6 


anlaytic, "+"  GMM-MUSCL, GMM-ENO,  dx=dt=1/100 


(a)  Ax  —  At  —  1O0 


anlaytic, "+"  GMM-MUSCL, GMM-ENO,  dt=1/200 


(b)  Ax  —  100,  At  —  200 


anlaytic, "+"  GMM-MUSCL, GMM-ENO,  dt=1/300 


(c)  Ax  —  10(),  At  —  3f)0 

Fig  II.5  High  resolution  methods  MUSCL  and  ENO 
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